Relaxed Purifying Selection and Possibly High Rate of Adaptation in Primate Lineage-Specific Genes

Keywords: lineage-specific genes, nonsynonymous polymorphism, evolutionary age, adaptive selection

LS of Genes

LS describes how specifically a gene and orthologs of the gene are distributed on a given phylogeny (Cai, Woo, et al. 2006). If a gene and its orthologs are present in the species all belong to a single lineage, the gene is considered specific to this particular lineage. On the other hand, if a gene and its orthologs are present in all the species of all lineages on the phylogeny, the gene is a “common” gene not specific to any lineages. Most of genes, however, have a certain level of LS laying between those of two “extremes” scenario—they are present in some but not all species.

To calculate LS for human genes in regard to the primate lineage, we used the phylogeny of 11 eukaryotic species, including Homo sapiens, Pan troglodytes, Mus musculus, Bos Taurus, Gallus gallus, Xenopus tropicalis, Danio rerio, Ciona intestinalis, Drosophila melanogaster, Caenorhabditis elegans, and Saccharomyces cerevisiae (fig. 1A). The topology of the phylogeny, supported by a range of molecular and morphological data (Blair et al. 2002; Bourlat et al. 2006; Nikolaev et al. 2007; Dunn et al. 2008), was retrieved from the National Center for Biotechnology Information (NCBI) Taxonomy database (Wheeler et al. 2000). The distribution of each human gene (i.e., the presence–absence pattern of its orthologs) on the tree forms the phylogenetic profile of the gene. We obtained phylogenetic profiles of genes from PhyloPat v41 (Hulsen et al. 2006) (http://www.cmbi.ru.nl/pw/phylopat/). Phylopat uses information of orthologs predicted in Ensembl compara database (Birney et al. 2006) to construct phylogenetic profile for each gene according to the presence or absence of orthologs of the gene in other species. The software pipeline of Ensembl compara database collected gene pairs of the best reciprocal hits and best score ratio values from a WUBlastp or Smith–Waterman whole-genome comparisons and then created a graph of gene relations, followed by a clustering step. The clusters were then applied to build a multiple alignment and a phylogenetic tree, which is reconciled with the species tree. From each reconciled gene tree, the orthologous relationships were inferred (for details, see [Vilella et al. 2009] and http://www.ensembl.org/info/docs/compara/homology_method.html). In figure 1A, we represented the phylogenetic profiles of human genes with a string of 11 symbols: and indicate the presence and absence of ortholog in the corresponding species, respectively. In our data analysis, we only included genes whose string representation of phylogenetic profile belongs to one of ten “regular” given patterns, in which and are constitutively arranged, so that LS level can be unambiguously assigned to these regular patterns and hence to genes whose phylogenetic profiles match these patterns. We discarded 10,032 (out of 20,150) genes that showed “irregular” phylogenetic profiles.

PL of Genes

Phylostratum is a set of genes from an organism that coalesce to founder genes having common phylogenetic origin (Domazet-Loso et al. 2007). Using a phylostratigraphic approach, Domazet-Loso and Tautz (2008) assigned all human protein-coding genes into 19 phylostrata. Here, we describe the procedure they used to determine gene’s PL. Briefly, BlastP algorithm (E value cutoff 0.001) was used to search human proteins against the NCBI nonredundant (NR) database to determine the presence/absence of homologs. The 0.001 cutoff value presents a good compromise between specificity and sensitivity of sequence similarity search (Domazet-Loso and Tautz 2003). Before the sequence search, the NR database was cleaned up with respect to sequences of uncertain or missing taxonomic status, unreliable phylogenetic positions, and filled up with complete genomes that were absent in the database but otherwise were publicly available. In addition, TBlastN search was conducted against substantial trace and expressed sequence tags archives of Porifera, Cyclostomata, and Chondrichtyes as for these internodes complete annotated genomes were lacking.

We obtained PL of all human genes from the supplementary data of Domazet-Loso and Tautz (2008). The numbers of genes in different PL groups vary (supplementary fig. S1, Supplementary Material online). To facilitate data analysis, we merged 19 PL groups into nine bins (fig. 2A). Empirically, neighboring phylostrata (e.g., phylostrata 9–11) with few genes were merged into one bin. Phylostrata with more genes (e.g., phylostrata 8 and 6) were assigned as two separate bins. The “empirically populated binning” procedure ensured the difference of gene number between bins was less substantial (supplementary fig. S1, Supplementary Material online). We also binned genes using an alternative procedure—“equally populated binning.” We added a random variable, ϵ ∼ norm(0, 0.001), to PL value of each gene, making PL a continuous variable. The value of ϵ was small such as to not change the original rank of PL of genes substantially, but, by adding an ϵ, each gene obtained a distinct rank. To generate equally populated bins, we adjusted the width of bins so that the same number of genes would fall into each bin. Two binning procedures produce similar results in subsequent data analyses. We only report results obtained using the first procedure. Note that this analysis was performed with genes in Ensembl database v45, for which the PL values were originally estimated (Domazet-Loso and Tautz 2008).

Gene Sets and Divergence Rate

To estimate the divergence rate of human protein-coding genes, we retrieved coding sequences of 20,150 human genes and their corresponding orthologs in chimpanzee genome from Ensembl database (Birney et al. 2006; Hubbard et al. 2007). Nonsynonymous substitution rate (K_a) and synonymous substitution rate (K_S) for human–chimpanzee orthologous pairs were calculated using the counting method of Nei and Gojobori (1986) implemented in MBEToolbox (Cai et al. 2005; Cai, Smith, et al. 2006). Before comparing median values of K_a and K_s between different groups of genes, we removed data points of 70 ribosomal genes (that are extremely slow-evolving genes that lack divergence information), 228 sex-chromosomal genes (that are under different strength of selection compared with autosomal genes), and 1,997 pseudogenes listed in http://Pseudogene.org/ (Karro et al. 2007). We further removed data points with K_a ≥ 0.05 and/or K_S ≥ 0.05 from 523 human–chimpanzee gene pairs to reduce the problem introduced by misalignment of coding sequences. Similar overall results were obtained when K_a and K_S (or denoted as dN and dS) were calculated using the maximum-likelihood method (Goldman and Yang 1994) implemented in PAML (Yang 1997). Given the fact that chimpanzees are the most closely related species to humans and divergence between human–chimp ortholog pairs may be extremely low, it is possible that the low values of K_a and K_S for some genes may just be an artifact of the choice of species to do the comparison. To eliminate the concern, we recomputed K_a and K_S with human–macaque ortholog pairs. The sequences of the corresponding macaque orthologs were again obtained from Ensembl database (Birney et al. 2006; Hubbard et al. 2007). We obtained similar result as above (supplementary fig. S2, Supplementary Material online). In addition, instead of computing K_a and K_S for individual genes, we summed divergence sites of genes in the same age classes and calculated pooled K_a and K_S for each age class (same procedure as described below). We obtained qualitatively unchanged results (data not shown).

We also obtained the numbers of nonsynonymous and synonymous sites (N and S) and the numbers of nonsynonymous (D_n) and synonymous (D_s) differences in coding sequences of human–chimpanzee genes resequenced in the study of Bustamante et al. (2005), who applied exon-specific polymerase chain reaction amplification to 20,362 loci in 39 humans and one chimpanzee to obtain sequence variants in these regions. We summed up D_n (and D_s) and N (and S) for genes in the same group to calculate pooled K_a and K_S (i.e., K_a and K_S for groups, in which sequences of genes are essentially concatenated). To obtain the 95% confidence interval (CI) of pooled K_a and K_S, we used the bootstrapping approach. For each group of genes, we constructed 10,000 resamples of the observed gene sets (and of equal size to the observed gene sets), each of which was obtained by randomly sampling with replacement from the original genes in the group. The CIs of pooled K_a and K_S were then obtained using percentile method from values of pooled K_a and K_S for resamples.

SNPs and Polymorphic Consequences

We computed A* (the number of nonsynonymous single-nucleotide polymorphism [nSNP] per nonsynonymous site) as the ratio between total nSNPs and total number of nonsynonymous sites of genes in the same class, and S* (the number of synonymous SNPs [sSNPs] per synonymous site) as the ratio between total sSNPs and total number of synonymous sites of genes in the same class. To determine the nonsynonymous and synonymous status of SNPs, we mapped SNPs (from dbSNP or Perlegen) onto the coding regions of the longest transcripts of human genes using functions in PGEToolbox (Cai 2008). For Applera SNPs, we used the nonsynonymous and synonymous classification from the original study (Bustamante et al. 2005). We systematically repeated SNP-related analyses using three SNP data sets: 1) Applera SNPs in the 11,624 resequenced genes from the study of Bustamante et al. (2005), 2) SNPs in dbSNP build 126, and 3) Perlegen SNPs (Hinds et al. 2005). Applera SNPs were identified by resequencing and thus suffer from less severe ascertainment biases. For dbSNP data, we removed SNPs that are not validated (dbSNP category 0) and retained NCBI designated “double hit” or “submitter validated” SNPs, which are more likely to beveal. For Perlegen data, we used all SNPs ascertained in all populations. We obtained allele frequencies of Applera SNPs and Perlegen SNPs to perform allele frequency spectrum analysis.

Robustness Tests for SNP Analyses

To conduct robustness tests against several confounding factors, we identified singleton genes, masked CpG-related SNPs, grouped genes according to functional categories, and classified genes according to their local genomic guanine-cytosine (GC) content. To identify singleton genes, we retrieved the data set of predicted human paralogous genes from Ensembl database (Birney et al. 2006; Hubbard et al. 2007). Singleton genes are those that do not have any paralogs.

Single-nucleotide mutations at CpG sites are much more frequent than at other sites (Cooper and Krawczak 1989; Sved and Bird 1990). To control for the effect of CpG-related SNPs, we generated the subset of SNPs excluding all CpG-related SNPs. To do so, we downloaded neighboring nucleotides of each SNP from human genome sequence from Ensembl v41 (Birney et al. 2006). SNPs were considered potentially CpG related in cases when: 1) A/G SNPs proceeded by a C or 2) C/T SNPs followed by G (Webster and Smith 2004).

The association between Gene Ontology (GO) terms and individual genes was obtained from FatiGO (Al-Shahrour et al. 2004). We used genes whose functional annotation contains the same GO term and tested whether the pattern between K_a/K_S (or A*/S*) and LS remains the same as that obtained from all genes. We obtained the result for GO terms: “cellular physiological process,” “metabolism,” “regulation of cellular process,” “regulation of physiological process,” and “cell communication,” which were associated with at least 25 genes in all of LS groups.

To test the effect of regional nucleotide composition of genes, we obtained the isochore map of human genome from the supplementary table of Costantini et al. (2006). Genes within isochore family L1 and L2 were considered in the GC-poor regions; genes within isochore families H1, H2, and H3 were considered in the GC-rich regions.

We obtained global mRNA expression data from Gene Expression Atlas (http://wombat.gnf.org) (Su et al. 2004). We included normal adult samples from 54 NR tissue types in the analysis. The expression level of each probe set in a given tissue was calculated as the mean of log (base 2) signal intensities of all samples after GC-robust multi-array average normalization (Wu et al. 2004). When multiple probes were mapped onto the same gene, the probe with the highest expression level was used as the report probe for this gene. We calculated the mean expression level of a gene (aveExp) as the mean of log signal intensities of probe sets across all tissues. We also calculated the maximum expression level among all tissues (maxExp) and the heterogeneity of expression level across all tissues (hetExp) for all genes with expression data available according to our previous study (Cai et al. 2009).

Estimation of the Ratio of the Rate of Adaptive Nonsynonymous Substitution to Synonymous Substitution (ω_A)

The distribution of fitness effects (DFEs) of new deleterious mutations and ω_A were estimated by using the method of Eyre-Walker and Keightley (2009) implemented in the program DoFE v2.1 (http://www.lifesci.susx.ac.uk/home/Adam_Eyre-Walker/Website/Software.html). To make the input file, we compiled Applera data obtained from Bustamante et al. (2005) and Lohmueller et al. (2008), which contain the numbers of nucleotides divergent between human and chimp and the site frequency spectra (SFS) for sSNPs and nSNPs, respectively. Genes of different LS groups were separated into different sets. Each set comprises the numbers of selected and neutral divergence polymorphism sites, as well as SFS vectors. DoFE provides an option for excluding singletons in SFS. Our result was unaffected qualitatively using the analysis with this option, thus only the result obtained without excluding singletons is reported below.

Estimation of Selection Coefficient (γ) using MKPRF

We used the program MKPRF (http://cbsuapps.tc.cornell.edu/mkprf.aspx) to estimate the selection coefficient γ (=2N_es) of genes. The program MKPRF samples from the posterior distribution of parameters in the MKPRF models of (Bustamante et al. 2002) and (Barrier et al. 2003) using a Markov Chain Monte Carlo algorithm based on Poisson random field (PRF) theory (Sawyer and Hartl 1992). We used a subset of genes with at least two variable nonsynonymous sites in the alignment (i.e., P_n + D_n ≥ 2). Using exactly the same approach as implemented in Bustamante et al. (2005), we applied the nonhierarchical model by specifying the flag FIXED_VARIANCE = 1, such that a single Gaussian prior of γ with a mean of 0 and the standard deviation (SD) (σ) of 8 is set for all loci (Bustamante CD, personal communication). Slightly deleterious SNPs can lead to an underestimate of the rate of adaptive evolution because they contribute to polymorphism but rarely become fixed. The effects of these slightly deleterious mutations can be partially controlled by removing mutations segregating at low frequencies (Fay et al. 2001). To circumvent the problem of unequal levels of slightly deleterious polymorphisms present in genes of different LS groups, we used two related procedures: 1) removing all SNPs whose derived alleles are at low frequencies (derived allele frequency [DAF] < 0.15) and 2) subsampling nSNPs at low frequencies (DAF < 0.15), such that SFS distributions across all LS classes match each other. Before subsampling, we calculated the fractions of these low-frequency nSNPs (lfnSNPs) in all nSNPs for genes in ten LS groups. The nSNPs in LS 10 group genes have the lowest fraction of lfnSNPs. We thus kept lfnSNPs in LS 10 group unchanged and subsampled lfnSNPs in the rest of the LS groups. For each LS group 1–9, we computed how many lfnSNPs (ϕ%) should be removed in order to make the final fraction of lfnSNPs equal to that of LS 10. Then we randomly purged ϕ% of lfnSNPs from each of LS groups 1–9. As a result, allele frequency spectra of the subsampled nSNPs in all ten LS groups became similar to each other (supplementary fig. S3, Supplementary Material online). We also used two different DAF cutoffs (<0.05 and <0.20) to define lfnSNPs. The effect of subsampling on final MKPRF results does not change depending on the exact values of the DAF cutoffs.

Li et al. (2008) found that results of MKPRF are sensitive to the model, and the value of σ used to estimate γ values. To get a sense of the robustness of γ estimation, we reran our analyses using hierarchical model (FIXED_VARIANCE = 0) with all default parameters and also using nonhierarchical model with the values of σ set at 1, 4, and 16.

Classifications of Human–Chimp Genes

We first classified human genes into ten groups based on their LS (Cai, Woo, et al. 2006). Here, we considered the primate lineage containing human and chimpanzee in the phylogeny of total 11 eukaryotic species (fig. 1A). We used LS to quantify the extent to which orthologs of a human gene are distributed in the species close to the primate lineage. Human genes whose orthologs are present in few species closely related to human and chimp are more primate lineage-specific than those genes whose orthologs are present in more distantly distributed species. We only assessed genes of ten LS categories that have nonpatchy distributions (fig. 1A). The numbers of genes in the LS categories 1–10 are 737, 848, 355, 458, 1,377, 294, 361, 834, 219, and 1,817, respectively.

We then classified human genes based on their PL, which quantifies evolutionary age of genes in terms of the most diverged lineage in which the homologs of those genes can be detected using Blast (see Materials and Methods and [Domazet-Loso et al. 2007]). As noted in Domazet-Loso and Tautz (2008), genes that share a particular protein domain will have the same PL based on when this domain emerged first, even though a particular gene may have evolved later (e.g., due to gene duplication or exon shuffling). Because a protein domain is usually linked to a certain function, PL is used to trace the origin of this function, irrespectively of the further origin of paralogs. We classified human genes into nine PL groups (Materials and Methods) (fig. 2A). LS and PL measures two proprieties of phylogenetic distribution of genes—breadth and depth, which are related to each other. Lineage-specific genes with high LS levels necessarily have high PL levels and vice versa. Indeed, the Spearman’s correlation coefficient between LS and PL is 0.65, and it is highly significant (P << 0.001).

Finally, we identified genes that were present in human, chimp, and yeast but that varied in their presence and absence in the “intermediate” lineages between human/chimp and yeast (fig. 3A). We classified these genes into four groups according to the number of GLs: (0) genes that have never been lost, (1) genes that have been lost once, (2) twice, and (3) three times or more (fig. 3A). These genes belong to the same PL group but vary in the propensity for loss. In such a setting, all human/chimp genes are successfully detected in the distantly related species (i.e., yeast), so it is much less likely that the patterns of presence or absence of genes among less divergent intermediate species are due to the failure of detecting genes by Blast.

Note that in every case, we retained only those human genes that have clear orthologs present in the chimpanzee genome. In this way, for all genes in the study, the rates of molecular evolution are estimated through human–chimp comparison in a consistent and reliable manner.

Higher Rates of Protein Evolution in Lineage-Specific Genes

We use the nonsynonymous substitution rate (K_a) and the ratio of nonsynonymous to synonymous substitution rate (K_a/K_S) between the human and chimpanzee orthologs to measure the rate of protein evolution. K_a and K_S values for the genes in each group (LS, PL, or GL) were calculated in two ways. First, we calculated K_a and K_S values for the individual genes and estimated the median value within each group. Second, we concatenated the sequences of all genes in the same group and calculated pooled K_a and K_S values for each group in this manner. Two approaches gave essentially identical values. For brevity, we report the results of the pooled K_a and K_S, which are based on the divergence data from the study of Bustamante et al. (2005), unless stated otherwise.

The K_a and K_a/K_S values vary significantly among different LS, PL, and GL groups (Kruskal–Wallis [KW] test, P ≤ 0.001 in all cases) (figs. 1B, 2B, and 3B). K_a and K_a/K_S are positively correlated with the LS and PL and negatively with GL values (Spearman’s r = 0.503 and 0.507 [LS] and r = 0.215 and 0.206 [PL] for K_a and K_S, respectively, P < 0.001), whereas Ks shows either no correlation (LS and GL, Spearman’s P > 0.05) or significant but very weak correlation (PL, Spearman’s ρ = 0.064, P = 1.19 × 10⁻¹³). Given the fact that chimpanzees are the most closely related species to humans and divergence between human–chimp ortholog pairs may be extremely low, it is possible that the low values of K_a and K_S for some genes may just be an artifact of the choice of species to do the comparison. To eliminate the concern, we recomputed K_a and K_S with human–macaque ortholog pairs. The sequences of the corresponding macaque orthologs were again obtained from Ensembl database (Birney et al. 2006; Hubbard et al. 2007). We obtained similar result as above (supplementary fig. S2, Supplementary Material online). These results are consistent with previous reports in humans or species in other domains (Krylov et al. 2003; Alba and Castresana 2005; Cai, Woo, et al. 2006).

Robustness to Potential Inability to Detect Fast-Evolving Genes in Distant Lineages

Inability of Blast to detect orthologs of fast-evolving genes in distant lineages can in principle contribute to the inverse relationship between the protein-divergence rate and LS (or PL) (Elhaik et al. 2006). This problem is unlikely to be severe, as simulations of the evolution of protein sequences with the same rates and among-site heterogeneity as those estimated from real mammalian protein-coding genes demonstrated that most functional genes could be detected by Blast in comparisons of even very distantly related genomes (e.g., fungi or plants vs. mammals) (Alba and Castresana 2007). The simulated procedure is exactly the same as the one used for the determination of PL values and should have as much sensitivity as the one employed in the LS statistic. Furthermore, the GL statistic, which measures the rate of GL, should be the least susceptible to this problem. This is because in the case of GL, proteins are first detected between the most distant lineages (humans and yeast) and thus are unlikely to be missed in the comparisons of more closely related species. Finally, in our analyses, the strongest signal comes from the comparisons of fairly young genes (e.g., high LS groups, fig. 1) and thus should be more robust to this problem.

We also conducted an additional test. We divided all genes into two groups: the slower evolving group (K_a ≤ 0.007) and the faster evolving half (K_a > 0.007). The cutoff K_a = 0.007 is the 72th percentile of ordered K_a values, which was chosen to split genes into two groups in such a way that each group contains enough (>30) genes in each of 10 LS or 9 PL classes. The positive correlation between K_a/K_S and LS and PL values is evident in both data sets (P < 0.001 in both cases, fig. 4). This demonstrates that the detected pattern is not due to the unusual behavior of the fastest evolving genes, which are the likeliest genes to go undetected in distant comparisons.

Overall, we believe it is highly unlikely that the problem of detecting proteins in distantly related lineages is sufficiently severe to invalidate the described analyses. Note that the slow rate of protein evolution is evident in the older and more broadly distributed genes independently of whether we classify genes using the LS, PL, or GL statistics. All the analyses described below can be reproduced using any of the statistics. In the rest of the manuscript, for brevity, we only report the results derived using the LS classification.

Robustness to the Variation in the Levels of Gene Expression

The level of gene expression is strongly anticorrelated with the protein divergence rate (Drummond et al. 2006; Pal et al. 2006). It is therefore important to determine whether the correlation between K_a or K_a/K_S and LS categories is due to the variation in expression levels of genes of different LS categories.

We consider three expression-related measures across all tissues: average expression (aveExp), the maximum expression (maxExp), and the heterogeneity of expression (hetExp) (see Materials and Methods). As expected, aveExp is anticorrelated with K_a/K_S (Spearman’s ρ = −0.256, P << 0.001) and K_a (Spearman’s ρ = −0.241, P << 0.001). The genes in higher LS categories have lower levels of expression (Spearman’s ρ = −0.230, P << 0.001). After controlling for expression level, the correlations between LS and K_a/K_S (or K_a) remain positive and highly significant [corr(K_a/K_S,LS|aveExp) = 0.238 and corr(K_a,LS|aveExp) = 0.222, both P << 0.001, Spearman partial correlation test]. Thus, correlation between LS and the rate of protein evolution is not entirely driven by lower expression levels of high LS genes. The other two variables, maxExp and hetExp, have similar relationships with LS and K_a/K_S as aveExp does, which is not unexpected given the strong correlations between both maxExp and hetExp with aveExp (Spearman’s ρ = 0.772 and 0.294, respectively [both P << 0.001]).

We conducted a linear multiple regression in the forward stepwise manner to examine the contributions of LS, aveExp, maxExp, and hetExp on the variation in log(K_a). The regression model defines log(K_a) as a function of all corresponding variables (X_LS, X_aveE, X_maxE, and X_hetE):

Table 1 gives the result of the modeling procedure. The final model gives a global R² of 0.131 (P < 0.001), that is, more than 13% of the variation in log(K_a) is explained by this model. During the construction of the final model, two predictors most highly correlated with log(K_a)—LS and aveExp—were kept. The remaining variables, which have minor roles in overall regression, were excluded. The standardized coefficients were examined to determine the relative importance of the significant predictors. LS contributes more to the model than does aveExp, as shown by its larger absolute standardized coefficient 0.350 and t statistic of 20.712, compared with values of 0.041 and 2.414 for aveExp. This analysis suggests that LS is the most relevant predictor of the rate of protein divergence.

Weaker Purifying Selection in Lineage-Specific Genes

The slower protein evolution of older and more broadly distributed genes is most likely due to stronger purifying selection acting on these genes. Indications for this should include a lower density of nSNPs and lower frequencies of derived alleles of these nSNPs. We use human SNP data and human/chimp divergence to investigate this prediction.

We analyze three SNP data sets: 1) SNPs in genes resequenced by Applera (Bustamante et al. 2005), 2) SNPs in NCBI dbSNP build 126, and 3) Perlegen SNPs (Hinds et al. 2005). Coding SNPs were split into sSNPs and nSNPs. For each LS category, we computed A* (the number of nSNP per nonsynonymous site) and S* (the number of sSNP per synonymous site) (see Materials and Methods). Figure 5 shows the results derived from all three data sets. None of the conclusions related to the assessment of the strength of purifying selection (measured by using the ratio A*/S*) change qualitatively depending on the investigated SNP data sets. Below we only describe Applera results.

First, we test whether older genes are subject to stronger purifying selection in humans. Values of A* and A*/S* correlate strongly and positively with the value of LS (Spearman’s ρ = 0.881 and 0.874, P < 0.005 and 0.001, respectively). The values of A*/S* increase almost 2-fold (0.31–0.57) from the lowest to the highest values of LS (fig. 5A). In addition, DAFs of SNPs in the genes from the lower LS groups are skewed toward rare alleles relative to that of higher LS genes (fig. 6). Specifically, the proportion of rare nSNPs (DAF < 0.15) is negatively correlated with LS values (Spearman’s ρ = −0.794, P < 0.01), whereas, in contrast, the proportion of rare sSNPs are not correlated with the LS values (P = 0.45, fig. 6). The proportion of rare nSNPs is significantly greater than the proportion of rare sSNPs for the LS groups 1 through 5 (P < 0.01 for all tests, G-test with Bonferroni correction) but not for the LS groups 6 through 10 (P > 0.05 for all tests, G-test). We also tried different DAF cutoffs (0.05, 0.1, and 0.2) as well as investigated the pattern derived from minor allele frequencies (MAFs). Regardless of the different cutoffs of DAF and the use of derived or MAF, results are similar (supplementary figs. S4 and S5, Supplementary Material online). Thus, compared with younger genes, older genes contain fewer nSNPs per site, and the frequencies of derived or minor alleles of these SNPs are lower. These results indicate that younger genes are subject to weaker purifying selection at the protein-coding level.

Robustness Tests for Polymorphism Patterns

We conducted several tests to demonstrate that the relationship between the proportion of A*/S* and LS cannot be explained entirely by a number of confounding factors. Specifically, we demonstrated that stronger purifying selection acting on older genes can be detected within subsets of the data defined by 1) whether a gene belongs to a particular functional (GO) gene group (supplementary table S1, Supplementary Material online), 2) whether a gene has a duplicate elsewhere in the human genome (supplementary table S2, Supplementary Material online), 3) whether SNPs are due to C to T transitions at CpG sites (supplementary table S3, Supplementary Material online), and 4) whether a gene resides in a genomic region of high or low GC content (supplementary table S4, Supplementary Material online) (for details, see Materials and Methods). The correlation between A*/S* and LS is unlikely due to different LS groups containing peculiar assemblages of genes defined by function, presence of a duplicate, number of CpG sites, or GC content.

We also verified that variation in levels of purifying selection in genes of different LS categories is not due to the variation in the level of gene expression, although A*/S* is strongly influenced by the expression pattern and breadth (Osada 2007). The partial correlation between A*/S* and LS is positive and significant after controlling for the different average expression levels of genes in different LS categories (i.e., corr(A*/S*,LS|aveExp) = 0.830, P < 0.001, partial Spearman correlation test with Applera data).

Rate of Adaptive Evolution in Lineage-Specific Genes

The restricted phylogenetic distribution of young genes implies that many of them are probably involved in lineage-specific adaptive processes. So, in addition to experiencing weaker purifying selection, young genes might be expected to experience higher rates of positive selection. To test this conjecture, we first used the method of Eyre-Walker and Keightley (2009) implemented in DoFE v2.1 to estimate ω_A for each LS categories. The method of Eyre-Walker and Keightley (2009) attempts to correct for to the problem existing in previous methods (e.g., those of Fay et al. [2001]; Smith and Eyre-Walker [2002]; Welch [2006]) that may give downwardly biased estimation if there are slightly deleterious mutations that inflate polymorphism relative to divergence (Crow and Kimura 1970; McDonald and Kreitman 1991; Eyre-Walker and Keightley 1999; Eyre-Walker 2002; Eyre-Walker et al. 2002). This method also estimates the DFEs of new deleterious mutations from the polymorphisms data and then uses the inferred DFE to predict the numbers of substitutions originating from neutral and slightly deleterious mutations between two species.

Using Applera polymorphism data and average allele frequencies across all African American and European American individuals, we estimated values of ω_A for ten LS categories, 0.0397, 0.0756, −0.0485, 0.0093, 0.0833, −0.0471, −0.031, 0.0155, −0.0756, and 0.1372. Because ω_A values and the LS values do not correlated with each other (P = 0.865, Spearman correlation test), these results provide no evidence that young genes experience a higher rate of adaptive nucleotide substitutions.

The DFEs of new neutral and deleterious mutations are simultaneously inferred by DoFE. We compared the estimated fractions of mutations in different N_es ranges among LS categories (fig. 7). There is a wide range of fractions of amino acid-changing mutations that behave as effectively neutral (0 < N_es < 1) among LS categories, ranging from 15.7% for LS 2 to 52.3% for LS 9. The fraction of effectively neutral mutations is correlated with LS significantly (Kendall’s τ = 0.73, P < 0.005), that is, the younger the genes the more effectively neutral mutations they have proportionally. The combined fractions of mutations with N_es ranges 1–10 and 10–100 are for slightly deleterious mutations. The fraction is negatively correlated with LS (r = −0.76, P = 0.017, Pearson correlation test after excluding outlier data point at LS 5) (fig. 7 and supplementary fig. S6, Supplementary Material online). This is consistent with previous results that old genes have proportionally more slightly deleterious mutations.

We also use the MKPRF analysis (Sawyer and Hartl 1992; Bustamante et al. 2002, 2005) to assess the proportion of genes showing evidence of positive selection in each LS categories. First, we ran MKPRF with a nonhierarchical model and a SD (σ = 8) of Gaussian prior, using exactly the same settings as in Bustamante et al. (2005). For each gene, we estimated the value of population-effective selection coefficient, γ (=2N_es, where N_e is the effective population size and s is the selection coefficient in a Wright–Fisher genic selection model). The values of γ vary significantly among genes of different LS groups (P < 0.001, KW test) and are positively correlated with LS values (Spearman’s ρ = 0.706, P < 0.001). We also obtained the 95% CIs of γ. If a gene has its 95% CIs of γ completely above 0, the gene appears to have been evolving under positive selection. On the other hand, if the 95% CIs of γ are completely below 0, the gene appears to be under negative selection and have a high proportion of weakly deleterious nonsynonymous polymorphisms. We found that the proportion of positively selected genes (f_γ>0) increases with the increment of LS values (Spearman’s ρ = 0.927, P < 0.001; Kendall’s τ = 0.778, P < 0.001) (fig. 8A, left panel), and the proportion of negatively selected genes (f_γ<0) decreases with the increment of LS values (Spearman’s ρ = −0.924, P < 0.001; Kendall’s τ = −0.778, P < 0.001) (fig. 8A, right panel). We also calculated the ratios of the numbers of positively selected and negatively selected genes to the numbers of neutrally evolving genes (f-γ>0 and f-γ<0, respectively). Similar to f_γ>0 and f_γ<0, f-γ>0 correlates positively with the LS values (Spearman’s ρ = 0.927; P < 0.001; Kendall’s τ = −0.778, P < 0.001) (fig. 8A, left panel), and f-γ<0 correlates negatively with the LS values (Spearman’s ρ = −0.915; Kendall’s τ = −0.778, P < 0.001) (fig. 8A, right panel). Note that, as in Bustamante et al. (2005), we focused our analysis only on the potentially informative loci that contain enough polymorphism and/or divergence events to have a chance of showing signals of selection. Specifically, we included informative loci with P_n + D_n ≥ 2 (Bustamante et al. 2005). We also carried out MKPRF analysis in two ways: either analyzing all genes together or analyzing genes from each LS group separately. The results remain virtually unchanged. Therefore, MKPRF analysis with full Applera SNPs suggested that younger genes experience a higher rate of adaptive evolution.

Slightly deleterious SNPs lead to an underestimation of the rate of adaptive evolution because they make a substantial contribution to polymorphism but fix at a much lower rate compared with neutral polymorphisms (Crow and Kimura 1970; Eyre-Walker and Keightley 1999; Eyre-Walker 2002; Eyre-Walker et al. 2002). From allele frequency analysis of nSNPs and the MKPRF analysis, we know that genes from lower LS groups have more slightly deleterious polymorphisms on average. This is indicated both by the higher proportion of rare nSNPs (fig. 6) and by the higher proportion of genes for which the MKPRF estimate of γ is negative for the low LS genes. Unlike the DoFE method of Eyre-Walker and Keightley (2009), in which the effect of slightly deleterious mutations are attempted to be controlled, MKPRF analysis per se does not control for this effect. Thus, the finding of a smaller proportion of genes experiencing positive selection in low LS group genes might be an artifact of the larger proportion of slightly deleterious polymorphisms in these genes.

We try to rule out this possibility using two approaches: 1) removing SNPs of low frequencies (Fay et al. 2001) prior to the analysis and 2) subsampling nSNPs to ensure that nSNPs in different LS groups have same frequency distribution and thus have the same bias (see Materials and Methods). The first procedure attempts to limit the effect of slightly deleterious SNPs by focusing on more frequent SNPs under the assumption that slightly deleterious SNPs are unlikely to reach intermediate frequencies (Fay et al. 2001). After removing the SNPs at less than 15% frequency and rerunning the MKPRF procedure, the positive correlation between LS and portion of genes under positive selection remains virtually unchanged (fig. 8B, Spearman corr(LS, f_γ>0) = 0.93, P < 0.0001, and corr(LS, f-γ>0) = 0.94, P < 0.0001). The second, subsampling procedure attempts to subsample the nSNPs such that the DFE of all of them is similar and is affected by the slightly deleterious SNPs both as little as possible, but more importantly, to the same extent across all the LS groups. We subsampled the nSNPs in all groups using the distribution in the LS 10 group. After subsampling, all LS groups have indistinguishably similar proportions of slightly deleterious polymorphisms (P ∼ 1, χ² test, supplementary fig. S1, Supplementary Material online). Note also that the proportions of rare SNPs are not different in the synonymous and nonsynonymous classes in the LS 10 group genes (P = 0.517, G-test) suggesting that the proportion of the slightly deleterious nSNPs in the LS 10 group genes is negligible. Therefore, subsampling should sharply reduce the influence of slightly deleterious nSNPs on the analysis overall. The proportion of genes containing a detectable number of slightly deleterious nSNPs no longer decreases with the increment of LS when MKPRF is carried out with the subsampled SNPs (fig. 8C). Importantly, the proportion of genes under positive selection remains higher in the higher LS groups (fig. 8C, Spearman corr(LS, f_γ>0) = 0.94, P < 0.0001, and corr(LS, f-γ>0) = 0.95, P < 0.0001).

Results of the MKPRF analysis might change depending on the prior and the model used in the analysis (Li et al. 2008). To explore these effects, we reran the MKPRF analysis using the nonhierarchical model with three additional σ values (1, 4, and 16). We found that although the absolute values of portion of genes in each age groups change, the positive correlation between LS and f_γ>0 (or f-γ>0) remains, especially in the cases of high σ values (e.g., σ = 16, Spearman corr(LS, f_γ>0) = 0.93, P = 0.0001, and corr(LS, f-γ>0) = 0.88, P = 0.008) (supplementary fig. S7, Supplementary Material online). We also ran MKPRF analysis using the hierarchical model (see Materials and Methods). The positive correlation between LS and f_γ>0 (or f-γ>0) is weaker in this case but remains significant after removing rare frequent SNPs (P = 0.004, Spearman correlation test) and after the subsampling (P = 0.003, Spearman correlation test) (supplementary fig. S8, Supplementary Material online). Overall, the MKPRF results do suggest that younger genes tend to experience positive selection more frequently, although to a modest degree. The faster evolution of younger genes appears to be attributable almost entirely to the weaker purifying selection acting on these genes.

1. Alba MM, Castresana J. Inverse relationship between evolutionary rate and age of mammalian genes. Mol Biol Evol. 2005;22:598–606. doi: 10.1093/molbev/msi045. [DOI] [PubMed] [Google Scholar]
2. Alba MM, Castresana J. On homology searches by protein Blast and the characterization of the age of genes. BMC Evol Biol. 2007;7:53. doi: 10.1186/1471-2148-7-53. [DOI] [PMC free article] [PubMed] [Google Scholar]
3. Al-Shahrour F, Diaz-Uriarte R, Dopazo J. FatiGO: a web tool for finding significant associations of Gene Ontology terms with groups of genes. Bioinformatics. 2004;20:578–580. doi: 10.1093/bioinformatics/btg455. [DOI] [PubMed] [Google Scholar]
4. Barrier M, Bustamante CD, Yu J, Purugganan MD. Selection on rapidly evolving proteins in the Arabidopsis genome. Genetics. 2003;163:723–733. doi: 10.1093/genetics/163.2.723. [DOI] [PMC free article] [PubMed] [Google Scholar]
5. Berglund AC, Sjolund E, Ostlund G, Sonnhammer EL. InParanoid 6: eukaryotic ortholog clusters with inparalogs. Nucleic Acids Res. 2008;36:D263–D266. doi: 10.1093/nar/gkm1020. [DOI] [PMC free article] [PubMed] [Google Scholar]
6. Birney E, et al. Ensembl 2006. Nucleic Acids Res. 2006;34:D556–D561. doi: 10.1093/nar/gkj133. [DOI] [PMC free article] [PubMed] [Google Scholar]
7. Blair JE, Ikeo K, Gojobori T, Hedges SB. The evolutionary position of nematodes. BMC Evol Biol. 2002;2:7. doi: 10.1186/1471-2148-2-7. [DOI] [PMC free article] [PubMed] [Google Scholar]
8. Bourlat SJ, et al. Deuterostome phylogeny reveals monophyletic chordates and the new phylum Xenoturbellida. Nature. 2006;444:85–88. doi: 10.1038/nature05241. [DOI] [PubMed] [Google Scholar]
9. Bustamante CD, et al. Natural selection on protein-coding genes in the human genome. Nature. 2005;437:1153–1157. doi: 10.1038/nature04240. [DOI] [PubMed] [Google Scholar]
10. Bustamante CD, et al. The cost of inbreeding in Arabidopsis. Nature. 2002;416:531–534. doi: 10.1038/416531a. [DOI] [PubMed] [Google Scholar]
11. Cai JJ. PGEToolbox: a Matlab toolbox for population genetics and evolution. J Hered. 2008;99:438–440. doi: 10.1093/jhered/esm127. [DOI] [PubMed] [Google Scholar]
12. Cai JJ, Borenstein E, Chen R, Petrov DA. Similarly strong purifying selection acts on human disease genes of all evolutionary ages. Genome Biol Evol. 2009;1:131–144. doi: 10.1093/gbe/evp013. [DOI] [PMC free article] [PubMed] [Google Scholar]
13. Cai JJ, Smith DK, Xia X, Yuen KY. MBEToolbox: a MATLAB toolbox for sequence data analysis in molecular biology and evolution. BMC Bioinformatics. 2005;6:64. doi: 10.1186/1471-2105-6-64. [DOI] [PMC free article] [PubMed] [Google Scholar]
14. Cai JJ, Smith DK, Xia X, Yuen KY. MBEToolbox 2: an enhanced MATLAB toolbox for molecular biology and evolution. Evol Bioinform. 2006;2:189–192. [PMC free article] [PubMed] [Google Scholar]
15. Cai JJ, Woo PC, Lau SK, Smith DK, Yuen KY. Accelerated evolutionary rate may be responsible for the emergence of lineage-specific genes in ascomycota. J Mol Evol. 2006;63:1–11. doi: 10.1007/s00239-004-0372-5. [DOI] [PubMed] [Google Scholar]
16. Cooper DN, Krawczak M. Cytosine methylation and the fate of CpG dinucleotides in vertebrate genomes. Hum Genet. 1989;83:181–188. doi: 10.1007/BF00286715. [DOI] [PubMed] [Google Scholar]
17. Costantini M, Clay O, Auletta F, Bernardi G. An isochore map of human chromosomes. Genome Res. 2006;16:536–541. doi: 10.1101/gr.4910606. [DOI] [PMC free article] [PubMed] [Google Scholar]
18. Crow JF, Kimura M. An introduction to population genetics theory. New York: Harper & Row; 1970. [Google Scholar]
19. Daubin V, Ochman H. Bacterial genomes as new gene homes: the genealogy of ORFans in E. coli. Genome Res. 2004;14:1036–1042. doi: 10.1101/gr.2231904. [DOI] [PMC free article] [PubMed] [Google Scholar]
20. Domazet-Loso T, Brajkovic J, Tautz D. A phylostratigraphy approach to uncover the genomic history of major adaptations in metazoan lineages. Trends Genet. 2007;23:533–539. doi: 10.1016/j.tig.2007.08.014. [DOI] [PubMed] [Google Scholar]
21. Domazet-Loso T, Tautz D. An evolutionary analysis of orphan genes in Drosophila. Genome Res. 2003;13:2213–2219. doi: 10.1101/gr.1311003. [DOI] [PMC free article] [PubMed] [Google Scholar]
22. Domazet-Loso T, Tautz D. An ancient evolutionary origin of genes associated with human genetic diseases. Mol Biol Evol. 2008;25:2699–2707. doi: 10.1093/molbev/msn214. [DOI] [PMC free article] [PubMed] [Google Scholar]
23. Drummond DA, Raval A, Wilke CO. A single determinant dominates the rate of yeast protein evolution. Mol Biol Evol. 2006;23:327–337. doi: 10.1093/molbev/msj038. [DOI] [PubMed] [Google Scholar]
24. Dunn CW, et al. Broad phylogenomic sampling improves resolution of the animal tree of life. Nature. 2008;452:745–749. doi: 10.1038/nature06614. [DOI] [PubMed] [Google Scholar]
25. Elhaik E, Sabath N, Graur D. The “inverse relationship between evolutionary rate and age of mammalian genes” is an artifact of increased genetic distance with rate of evolution and time of divergence. Mol Biol Evol. 2006;23:1–3. doi: 10.1093/molbev/msj006. [DOI] [PubMed] [Google Scholar]
26. Eyre-Walker A. Changing effective population size and the McDonald-Kreitman test. Genetics. 2002;162:2017–2024. doi: 10.1093/genetics/162.4.2017. [DOI] [PMC free article] [PubMed] [Google Scholar]
27. Eyre-Walker A, Keightley PD. High genomic deleterious mutation rates in hominids. Nature. 1999;397:344–347. doi: 10.1038/16915. [DOI] [PubMed] [Google Scholar]
28. Eyre-Walker A, Keightley PD. Estimating the rate of adaptive molecular evolution in the presence of slightly deleterious mutations and population size change. Mol Biol Evol. 2009;26:2097–2108. doi: 10.1093/molbev/msp119. [DOI] [PubMed] [Google Scholar]
29. Eyre-Walker A, Keightley PD, Smith NG, Gaffney D. Quantifying the slightly deleterious mutation model of molecular evolution. Mol Biol Evol. 2002;19:2142–2149. doi: 10.1093/oxfordjournals.molbev.a004039. [DOI] [PubMed] [Google Scholar]
30. Fay JC, Wyckoff GJ, Wu CI. Positive and negative selection on the human genome. Genetics. 2001;158:1227–1234. doi: 10.1093/genetics/158.3.1227. [DOI] [PMC free article] [PubMed] [Google Scholar]
31. Fisher RA. The genetical theory of natural selection. Oxford: The Clarendon Press; 1939. [Google Scholar]
32. Goldman N, Yang Z. A codon-based model of nucleotide substitution for protein-coding DNA sequences. Mol Biol Evol. 1994;11:725–736. doi: 10.1093/oxfordjournals.molbev.a040153. [DOI] [PubMed] [Google Scholar]
33. Heinen TJ, Staubach F, Haming D, Tautz D. Emergence of a new gene from an intergenic region. Curr Biol. 2009;19:1527–1531. doi: 10.1016/j.cub.2009.07.049. [DOI] [PubMed] [Google Scholar]
34. Hinds DA, et al. Whole-genome patterns of common DNA variation in three human populations. Science. 2005;307:1072–1079. doi: 10.1126/science.1105436. [DOI] [PubMed] [Google Scholar]
35. Hubbard TJ, et al. Ensembl 2007. Nucleic Acids Res. 2007;35:D610–D617. doi: 10.1093/nar/gkl996. [DOI] [PMC free article] [PubMed] [Google Scholar]
36. Hulsen T, de Vlieg J, Groenen PM. PhyloPat: phylogenetic pattern analysis of eukaryotic genes. BMC Bioinformatics. 2006;7:398. doi: 10.1186/1471-2105-7-398. [DOI] [PMC free article] [PubMed] [Google Scholar]
37. Johnson ME, et al. Positive selection of a gene family during the emergence of humans and African apes. Nature. 2001;413:514–519. doi: 10.1038/35097067. [DOI] [PubMed] [Google Scholar]
38. Karro JE, et al. Pseudogene.org: a comprehensive database and comparison platform for pseudogene annotation. Nucleic Acids Res. 2007;35:D55–D60. doi: 10.1093/nar/gkl851. [DOI] [PMC free article] [PubMed] [Google Scholar]
39. Kasuga T, Mannhaupt G, Glass NL. Relationship between phylogenetic distribution and genomic features in Neurospora crassa. PLoS One. 2009;4:e5286. doi: 10.1371/journal.pone.0005286. [DOI] [PMC free article] [PubMed] [Google Scholar]
40. Khalturin K, Hemmrich G, Fraune S, Augustin R, Bosch TC. More than just orphans: are taxonomically-restricted genes important in evolution? Trends Genet. 2009;25:404–413. doi: 10.1016/j.tig.2009.07.006. [DOI] [PubMed] [Google Scholar]
41. Kimura M. The neutral theory of molecular evolution. New York: Cambridge University Press, Cambridge [Cambridgeshire]; 1983. [Google Scholar]
42. Knowles DG, McLysaght A. Recent de novo origin of human protein-coding genes. Genome Res. 2009;19:1752–1759. doi: 10.1101/gr.095026.109. [DOI] [PMC free article] [PubMed] [Google Scholar]
43. Kouprina N, et al. The SPANX gene family of cancer/testis-specific antigens: rapid evolution and amplification in African great apes and hominids. Proc Natl Acad Sci U S A. 2004;101:3077–3082. doi: 10.1073/pnas.0308532100. [DOI] [PMC free article] [PubMed] [Google Scholar]
44. Krylov DM, Wolf YI, Rogozin IB, Koonin EV. Gene loss, protein sequence divergence, gene dispensability, expression level, and interactivity are correlated in eukaryotic evolution. Genome Res. 2003;13:2229–2235. doi: 10.1101/gr.1589103. [DOI] [PMC free article] [PubMed] [Google Scholar]
45. Kuo CH, Kissinger JC. Consistent and contrasting properties of lineage-specific genes in the apicomplexan parasites Plasmodium and Theileria. BMC Evol Biol. 2008;8:108. doi: 10.1186/1471-2148-8-108. [DOI] [PMC free article] [PubMed] [Google Scholar]
46. Li YF, Costello JC, Holloway AK, Hahn MW. “Reverse ecology” and the power of population genomics. Evolution. 2008;62:2984–2994. doi: 10.1111/j.1558-5646.2008.00486.x. [DOI] [PMC free article] [PubMed] [Google Scholar]
47. Lohmueller KE, et al. Proportionally more deleterious genetic variation in European than in African populations. Nature. 2008;451:994–997. doi: 10.1038/nature06611. [DOI] [PMC free article] [PubMed] [Google Scholar]
48. McDonald JH, Kreitman M. Adaptive protein evolution at the Adh locus in Drosophila. Nature. 1991;351:652–654. doi: 10.1038/351652a0. [DOI] [PubMed] [Google Scholar]
49. Nei M, Gojobori T. Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol Biol Evol. 1986;3:418–426. doi: 10.1093/oxfordjournals.molbev.a040410. [DOI] [PubMed] [Google Scholar]
50. Nikolaev S, et al. Early history of mammals is elucidated with the ENCODE multiple species sequencing data. PLoS Genet. 2007;3:e2. doi: 10.1371/journal.pgen.0030002. [DOI] [PMC free article] [PubMed] [Google Scholar]
51. Orr HA. The population genetics of adaptation: the adaptation of DNA sequences. Evolution. 2002;56:1317–1330. doi: 10.1111/j.0014-3820.2002.tb01446.x. [DOI] [PubMed] [Google Scholar]
52. Osada N. Inference of expression-dependent negative selection based on polymorphism and divergence in the human genome. Mol Biol Evol. 2007;24:1622–1626. doi: 10.1093/molbev/msm080. [DOI] [PubMed] [Google Scholar]
53. Pal C, Papp B, Lercher MJ. An integrated view of protein evolution. Nat Rev Genet. 2006;7:337–348. doi: 10.1038/nrg1838. [DOI] [PubMed] [Google Scholar]
54. Sawyer SA, Hartl DL. Population genetics of polymorphism and divergence. Genetics. 1992;132:1161–1176. doi: 10.1093/genetics/132.4.1161. [DOI] [PMC free article] [PubMed] [Google Scholar]
55. Smith NG, Eyre-Walker A. Adaptive protein evolution in Drosophila. Nature. 2002;415:1022–1024. doi: 10.1038/4151022a. [DOI] [PubMed] [Google Scholar]
56. Su AI, et al. A gene atlas of the mouse and human protein-encoding transcriptomes. Proc Natl Acad Sci U S A. 2004;101:6062–6067. doi: 10.1073/pnas.0400782101. [DOI] [PMC free article] [PubMed] [Google Scholar]
57. Sved J, Bird A. The expected equilibrium of the CpG dinucleotide in vertebrate genomes under a mutation model. Proc Natl Acad Sci U S A. 1990;87:4692–4696. doi: 10.1073/pnas.87.12.4692. [DOI] [PMC free article] [PubMed] [Google Scholar]
58. Toll-Riera M, et al. Origin of primate orphan genes: a comparative genomics approach. Mol Biol Evol. 2009;26:603–612. doi: 10.1093/molbev/msn281. [DOI] [PubMed] [Google Scholar]
59. Toll-Riera M, Castelo R, Bellora N, Alba MM. Evolution of primate orphan proteins. Biochem Soc Trans. 2009;37:778–782. doi: 10.1042/BST0370778. [DOI] [PubMed] [Google Scholar]
60. Toll-Riera M, Castresana J, Alba MM. Accelerated evolution of genes of recent origin. In: Pontarotii P, editor. Evolutionary biology from concept to application. Berlin (Germany): Springer; 2008. [Google Scholar]
61. Vilella AJ, et al. EnsemblCompara GeneTrees: complete, duplication-aware phylogenetic trees in vertebrates. Genome Res. 2009;19:327–335. doi: 10.1101/gr.073585.107. [DOI] [PMC free article] [PubMed] [Google Scholar]
62. Wang W, et al. Origin and evolution of new exons in rodents. Genome Res. 2005;15:1258–1264. doi: 10.1101/gr.3929705. [DOI] [PMC free article] [PubMed] [Google Scholar]
63. Webster MT, Smith NG. Fixation biases affecting human SNPs. Trends Genet. 2004;20:122–126. doi: 10.1016/j.tig.2004.01.005. [DOI] [PubMed] [Google Scholar]
64. Welch JJ. Estimating the genomewide rate of adaptive protein evolution in Drosophila. Genetics. 2006;173:821–837. doi: 10.1534/genetics.106.056911. [DOI] [PMC free article] [PubMed] [Google Scholar]
65. Wheeler DL, et al. Database resources of the National Center for Biotechnology Information. Nucleic Acids Res. 2000;28:10–14. doi: 10.1093/nar/28.1.10. [DOI] [PMC free article] [PubMed] [Google Scholar]
66. Wilson AC, Carlson SS, White TJ. Biochemical evolution. Annu Rev Biochem. 1977;46:573–639. doi: 10.1146/annurev.bi.46.070177.003041. [DOI] [PubMed] [Google Scholar]
67. Wolf YI, Novichkov PS, Karev GP, Koonin EV, Lipman DJ. The universal distribution of evolutionary rates of genes and distinct characteristics of eukaryotic genes of different apparent ages. Proc Natl Acad Sci U S A. 2009;106:7273–7280. doi: 10.1073/pnas.0901808106. [DOI] [PMC free article] [PubMed] [Google Scholar]
68. Wu Z, Irizarry RA, Gentleman R, Martinez-Murillo F, Spencer F. A model based background adjustment for oligonucleotide expression arrays. J Am Stat Assoc. 2004;99:909–917. [Google Scholar]
69. Yang Z. PAML: a program package for phylogenetic analysis by maximum likelihood. Comput Appl Biosci. 1997;13:555–556. doi: 10.1093/bioinformatics/13.5.555. [DOI] [PubMed] [Google Scholar]