Molecular evolution of genes in avian genomes

Pairwise comparison of the chicken and zebra finch protein-coding gene sets

We identified 11,225 1:1 orthologs from the pairwise comparison of all protein-coding genes in the chicken and zebra finch draft genome sequences. This corresponds to 60 to 65% of the total number of genes in the avian genome [10]. The overall degree of neutral divergence, as approximated by the rate of synonymous substitution (d _S ) from 1,000 random sets of 150 genes [25], between these two bird species was 0.418 (95% confidence interval = 0.387 to 0.458). The overall ω (d _N /d _S ) in the pairwise comparison was 0.152 (95% confidence interval = 0.127 to 0.179).

Lineage-specific rates of evolution

d _S was significantly (8%) higher in the zebra finch (0.213) than in the chicken lineage (0.197; P < 2.2 × 10^-16, Wilcoxon signed rank test; Table 1), indicating a difference in the molecular clock of these two parallel lineages. d _S of the ancestral bird lineage was higher (0.260) than in the two terminal branches, which is not unexpected given the estimated divergence times. The divergence at fourfold degenerate sites showed the same trend, and was highest in the ancestral bird lineage (mean of 1 Mb intervals = 0.239), and higher in zebra finch (0.199) than in chicken (0.172). We estimated lineage-specific mutation rates by dividing the divergence at fourfold degenerate sites with the estimated age of lineages according to the divergence times given above. We found that the mutation rate was lower in the ancestral bird lineage (1.23 × 10^-9 site^-1 year^-1)than in both the chicken lineage (1.91 × 10^-9 site^-1 year^-1; P < 2 × 10^-16) and the zebra finch lineage (2.21 × 10^-9 site^-1 year^-1; P < 2 × 10^-16), and that the rate in the chicken lineage was significantly lower than the rate in the zebra finch lineage (P < 1 × 10^-5).

The zebra finch lineage had a significantly higher overall ω than the chicken lineage (0.133 versus 0.121; P < 2.2 × 10^-16, Wilcoxon signed rank test). Just as for divergence, there was a strong correlation between individual ω values of 1:1 chicken and zebra finch orthologs (r² = 0.338, P < 2 × 10^-16). A corresponding analysis for 7,789 human and mouse orthologs (included in the 8,384 genes from multiple-species alignments) revealed a similarly strong correlation (r² = 0.359, P < 2 × 10^-16). Moreover, we also found a similar strength of correlation in gene-wise ω values estimated for orthologs from the bird lineage (chicken and zebra finch) with the mammalian lineage (human and mouse lineages; r² = 0.325, P < 2 × 10^-16). The gene-wise correlations between ω values for the ancestral bird lineage (which had an overall ω of 0.110) and chicken (r² = 0.178, P < 2 × 10^-16) and zebra finch (r² = 0.170, P < 2 × 10^-16), respectively, were weaker.

Adaptive evolution of genes in the avian genome

We next sought to identify genes, and the functional categories these genes are associated with, that are candidates for being involved with lineage-specific adaptations during avian evolution. We considered the ancestral bird lineage as well as the terminal chicken and zebra finch lineages separately, and posed three specific questions.

First, which genes have evolved most rapidly in avian lineages (high ω values), indicative of either adaptive evolution or relaxed selective constraint? For this question we used a likelihood ratio test to determine which genes had a significantly higher ω value than the mean of all genes in the genome. These genes are referred to as rapidly evolving bird (REB) genes. We used this approach rather than simply selecting, for example, the top 5% or 10% of genes sorted by ω value since the confidence in ω values is dependent of alignment length and the number of substitutions within a particular gene.

Second, which genes have evolved more rapidly in avian lineages than in other amniote lineages (mammals and lizard)? Here we used a branch model in PAML to determine which genes had a significantly higher ω in avian lineages than in other branches of the tree corresponding to our data. These genes are referred to as more rapidly evolving in birds (MREB).

Third, which genes show evidence of containing codons that have been subject to positive selection (referred to as PS genes) during avian evolution? For this third question we used a branch-site model in PAML to identify genes containing positively selected codons with ω higher than 1.

If positively selected codons are evenly distributed across genes and the power to detect such codons is more or less constant, then the likelihood of detecting genes containing positively selected codons will correlate with alignment length. Consistent with this, three out of three unique overrepresented GO terms from the list of positively selected genes in the ancestral bird branch have longer mean alignment length than genes with other GO terms (P < 0.001, Wilcoxon rank sum test). However, the overrepresented GO terms from the list of positively selected genes in the chicken lineage have actually shorter mean alignment length than genes with other GO terms, with marginal significance (P = 0.093). This warrants further investigation, from both methodological and biological points of view.

As a comparison, we tested for over-represented GO terms among positively selected mammalian genes and genes evolving significantly faster in mammals than in birds (Table S2 in Additional file 1). However, using the same criteria as applied to the lists of avian genes, no GO term was significantly over-represented in the mammalian lists.

Adaptive evolution of neurological genes

We then compared the number of genes classified as associated with 'learning', 'neurogeneration' or 'neurodegeneration' that were found to be positively selected in either the chicken or zebra finch lineage (that is, excluding genes that were positively selected in both lineages). Interestingly, for each OMIM term the number of unique positively selected genes was significantly higher in zebra finch than in chicken (Table 5; 10 versus 5, 27 versus 15, and 16 versus 8, respectively). This indicates that the songbird lineage has experienced more frequent adaptive evolution of genes relating to cognitive functions than the galliform lineage.

The relationship between selection and recombination

We sought to elucidate how the intensity of selection and/or the influence of genetic drift, manifested in ω, vary across the avian genome. The potential influence of recombination on ω was of particular interest since the rate of recombination is unusually heterogeneous within both the chicken [31] and zebra finch [32] genomes, and probably so for birds in general. Such heterogeneity could set the stage for recombination affecting the efficacy of selection and thereby ω, as predicted by evolutionary theory [33] but for which there is limited empirical support [34–38].

As a starting point for these analyses we first noted that there was a weak positive correlation between ω estimated for 1 Mb intervals and chromosome size in zebra finch (Figure 1; r² = 0.055, P = 6 × 10^-11) and chicken (r² = 0.029, P = 3 × 10^-6). This confirms similar observations made for a small set of chicken-turkey orthologs [11] as well as for chicken-human orthologs [10], although the effect we detected here with much larger data sets was considerably weaker than indicated by those previous studies. There was a strong negative correlation between the mean divergence of fourfold degenerate sites of 1 Mb intervals and chromosome size (Figure 2; r² = 0.153 in zebra finch and r² = 0.140 in chicken, P < 2 × 10^-16 in both cases). These correlations were not limited to the dichotomy of macrochromosomes versus microchromosomes (data not shown); indeed, for many birds chromosome size shows a relatively continuous distribution without a clear distinction between macrochromosomes and microchromosomes [7].

Over-represented gene ontologies among positively selected or rapidly evolving genes

With draft sequences now available for two avian genomes it is possible to study the role of natural selection in shaping individual gene sequences during avian evolution. An impetus for our study was thus to identify genes and gene categories that have been important for adaptive character evolution in a vertebrate lineage. Clearly, there are many morphological, physiological and behavioral phenotypes that distinguish birds and mammals. A comparative genomic approach has the potential to contribute towards the identification of the genetic basis of these differences [58].

Basic characteristics of birds such as feathers, flight and hollow bones evolved prior to the split of the chicken and zebra finch lineages. The genetic novelties underlying these phenotypes should thus have started to appear in an ancestral lineage. As discussed above, the ancestral bird branch in the phylogenetic tree formed by our data corresponds mostly to non-avian dinosaurs of the order Saurischia, suborder Theropoda. Genes or gene categories identified as positively selected or rapidly evolving in this branch may thus be related to phenotypic evolution in non-avian dinosaurs rather than in modern birds. On the other hand, many bird-like features may have started to emerge already for non-avian dinosaurs.

The two GO terms found to be over-represented among genes evolving under positive selection in the ancestral bird lineage, calcium ion binding and cell adhesion, largely represent an overlapping set of genes. Most of these genes (Table S1 in Additional file 1) encode transmembrane cadherins that play a critical role in cell-cell adhesion in tissue structures. One of these cadherins, protocadherin-15, is expressed in retina and we note that another positively selected calcium ion binding gene, Crumbs homolog 1, is involved with photoreceptor morphogenesis in retina; mutations in the human ortholog cause retinitis pigmentosa type 12 [59]. The visual ability of birds is superior to other vertebrates and the molecular adaptations underlying this phenotype are likely to have been driven by positive selection.

In the chicken lineage the term anion transmembrane transporter activity was over-represented among positively selected genes. The genes annotated with this term include solute carriers and ion channels involved with basic cell signaling processes, for example, in neurotransmission. In the zebra finch lineage the term microtubule cytoskeleton was over-represented among genes evolving faster in this lineage than in other branches of the amniote tree. The majority of these are kinesins and other genes involved with mitosis/meiosis, sperm motility, centrosome formation and synapse function.

It should be stressed that we inferred positive selection in lineages corresponding to nearly 100 million years or more of evolution and that large numbers of genes were uncovered by these analyses. This is likely to reduce the power of detecting enriched GO terms due to dilution and failure to capture temporal episodes of adaptive evolution. Moreover, given that our data were defined by a common set of 1:1 orthologous genes found in birds, a lizard and mammals, the analysis did not include lineage-specific genes that may be particularly responsive to positive selection. These aspects are probably of relevance to the somewhat surprising observation that no significantly over-represented GO terms were found among positively selected or rapidly evolving mammalian genes. This is seemingly at odds with previous work in primates that frequently have revealed categories such as sensory perception, immune defence, apoptosis and spermatogenesis to be enriched among positively selected genes [17, 18, 60–62]. In birds, there have recently been large-scale efforts toward transcriptome sequencing of several species, including songbirds [63]. These data will allow study of the molecular evolution of genes in much shorter branches of the avian phylogenetic tree than is currently possible with complete genome sequences, which is only available for chicken and zebra finch.

Zebra finch and positive selection in neurological genes

The zebra finch communicates through learned vocalizations ('songs'). Only the male zebra finch produces learned song, and he learns this song by copying an adult tutor during a critical period in juvenile development. Experimental work in zebra finch has demonstrated the localization and character of neural circuits involved in developmental song learning and adult singing [64–67], with dynamic regulation of brain gene expression in response to singing and song experience [68–76]. Fifty-eight genes with known roles in learning, neurogenesis or neurodegeneration in humans show evidence of positive selection in the zebra finch lineage. Of these, nine (15%; Table 6) are also expressed differentially in the song control system, either at higher or lower levels than in the surrounding brain tissue, according to the study of Lovell et al. [30]. In comparison, only 2% (390 out of 17,214) unique brain-derived cDNA probes on that microarray gave differential hybridization signals in the song control system. We note that five of the nine genes encode proteins involved in cell surface and synaptic signaling: voltage-dependent L-type calcium channel subunit alpha-1D (CACNA1D), G protein-coupled receptor 98 precursor (GPR98), glutamate receptor, ionotropic AMPA 2 (GRIA2), glutamate receptor, metabotropic 1 (GRM1), and protein tyrosine phosphatase receptor type F (PTPRF). GRIA12 is also one of the ion channel genes that are suppressed in response to song playbacks as reported in Warren et al. [16].

Four of the 58 genes show evidence of accelerated evolution in the primate lineage: ASPM, GRIN2a, DRD2, and LHX2. Two of these have apparent roles in neurogenesis and neuronal development (ASPM and LHX2). In particular, ASPM (abnormal spindle-like microcephaly-associated) has been a focus of speculation with respect to the dramatic evolution of brain size in humans. Homozygous mutations in ASPM are a cause of primary microcephaly [77] and the gene shows evidence of positive selection in both the human lineage [78–80] and the ancestral lineage of the apes [81]. Songbirds have also experienced a relative increase in brain size compared to other avian lineages [82], with the notable emergence of the large and highly plastic nuclei of the song control system. However, enthusiasm for ASPM as a key factor in primate brain evolution has been tempered by findings that mutations in ASPM are not correlated with cognitive ability [83, 84] and by alternative roles for ASPM that might place it under selection more broadly, such as a role in ciliary function [85].

The other two neurological genes that are also accelerated in primates may be considered to have neuromodulatory functions that can directly affect learning, memory and behavior. DRD2 encodes the D2 subtype of the dopamine receptor. GRIN2a encodes a subunit of the N-methyl-D-aspartate (NMDA) receptor, a subtype of ionotropic glutamate-gated ion channel that has well-established roles in learning and brain plasticity (reviewed in [86]). A survey of GRIN2a sequences across primate species revealed a specific correlation between ω and home range size, which is taken to be a proxy for spatial memory [87]. Spatial memory is well developed in the songbird (passerine) lineage and is especially evident in food-caching species [88], a behavior that depends on NMDA receptor function [89]. Zebra finches are not studied as a food caching species, but their nomadic lifestyle implies a highly sophisticated spatial sense [90]. NMDA receptors have also been implicated in song learning and song control system neurophysiology [91, 92]. The rich diversity of songbird species and their adaptations should provide unusual opportunities for correlating NMDA receptor sequence evolution with specific behavioral and neurophysiological variations.

The strength of selection during avian evolution

The overall strength of selection as manifested in the genome-wide ratio of non-synonymous to synonymous substitution rates (ω) was similar in the chicken and zebra finch lineages (0.12 to 0.13), as well as in the ancestral bird lineage (0.11). These ratios are about half that reported among hominids and more similar to what is seen in the murid and dog lineages [62]. This may be taken to suggest that the rate of adaptive evolution and/or the rate of accumulation of slightly deleterious mutations have been lower in birds than in primates. However, it is increasingly appreciated that point estimates of mean ω can be misleading. Mean ω decreases with branch length and needs to be seen in a time trajectory framework rather than as a fixed quantity [93–95]. The apparent lineage-specific differences between hominids on the one side and murids, dogs and birds on the other may thus simply be accounted for by branch length. Future research will be needed to explore how branch length is best accounted for when comparing mean ω for different lineages. This will be important when addressing whether life history variables, such as the effective population size (N_e), correlate with mean ω. For example, such a correlation might be expected if slightly deleterious mutations contribute significantly to protein evolution as postulated by a nearly neutral model [96], giving rise to a negative relationship between mean ω and N_e [97].

The relationship between natural selection and recombination in avian genomes

Selection acts in each generation on alleles embedded within a particular genomic background. Due to recombination, selection will, over time, be able to favor or disfavor alleles at individual loci without affecting the rest of the genome. This comes with a caveat that when two loci are genetically linked, selection at one locus will affect the efficiency of selection at the other: the loci are said to interfere with each other. Theory predicts that the strength of interference should be related to the amount of recombination between the loci; this is the so-called Hill-Robertson effect [33]. Theoretical predictions on the consequence of Hill-Robertson interference on coding sequence evolution depend on the fitness distribution of segregating variants at non-synonymous sites [98, 99]. If slightly deleterious mutations segregate frequently in the population, directional selection at one locus will increase the probability of fixation of such mutations at linked loci. If beneficial alleles are common in the population, the probability of fixation of those mutations will be reduced at linked loci. These two scenarios are associated with opposing predictions for the correlation between recombination rate and ω; in the former case a negative relationship is expected whereas in the latter case a positive relationship is expected.

The strongest support for Hill-Robertson interference comes from regions devoid of recombination. For example, ω is generally high in the non-recombining sex chromosome, that is, the Y chromosome in systems with male heterogamety and the W chromosome in systems with female heterogamety [100–102]. However, it has been surprisingly difficult to find genome-wide empirical support for Hill-Robertson interference, and data are currently limited to studies in Drosophila [34, 35, 103, 104] and a recent study of humans failed to demonstrate a correlation between recombination rate and ω [37].

It is possible that the power for detecting a relationship between recombination and ω could be higher in bird systems because the rate of recombination is highly heterogeneous, at least within the two avian genomes for which detailed information is currently available on regional recombination rate variation. Specifically, there is a clear negative relationship between chromosome size and recombination rate [10, 31] following from an obligate recombination event per chromosomal arm. In chicken, the average per-chromosome recombination rate ranges from 2 centiMorgans (cM)/Mb up to 10 cM/Mb [10]. Moreover, there is significant within-chromosome variation in the rate of recombination with a strong 'telomere effect'. This is most readily seen in zebra finch, with rates close to 10 cM/Mb in terminal regions of the larger (>100 Mb) chromosomes while the central parts have rates as low as 0.1 cM/Mb; the latter is not just a 'centromere effect' because these recombination deserts cover up to 75% of the larger chromosomes [32].

We do not find support for an increased efficiency of directional selection in regions of high recombination. If anything, the data go in the opposite direction since there was a weak negative, yet significant relationship between ω and recombination rate in both chicken and zebra finch (r² < 0.01, after controlling for GC content and the amount of coding sequence); this is the direction predicted from the hypothesis of an accumulation of slightly deleterious mutations in regions with low recombination rate. One obvious explanation for this weak relationship is that both slightly deleterious and beneficial variants are common and that their opposing effects in Hill-Robertson interference largely cancel each other out. However, in the absence of simulations under different distributions of the fitness consequences of segregation mutations this remains an argument that is difficult to examine.

Another explanation relates to the fact that recombination rate and ω are measured on very different time scales. Recombination is recorded from pedigree data and thus reflects the rate in contemporary populations. Lineage-specific ω represents substitutions that have accumulated during, in this case, 90 million years of avian evolution. If the recombination landscape has changed frequently during the course of this time period, this may have weakened the signal of potential recombination effects on the pattern of efficacy of selection across the genome. There is limited knowledge on the evolutionary consistency of regional recombination rate variation [105]. At a local scale, recombination hot-spots are ephemeral in the human genome with little or no evidence for hot-spots at orthologous positions in the chimpanzee genome [106–108]. As indicated above, recombination rates in birds are strongly associated with chromosome features, with highly elevated rates in microchromosomes and in telomeric regions of larger chromosomes. Given the high degree of karyotype stability in birds [7], this may suggest that the recombination landscape has also remained relatively stable. Indeed, we have found that recombination rates in 1-Mb windows of the chicken and zebra finch genomes to be significantly correlated [32]. Moreover, the strong correlation observed between base composition (GC content) and current recombination rates in both chicken [31] and zebra finch [32] is consistent with a conserved pattern of recombination rate variation, at least under the scenario that recombination drives the long-term evolution of base composition (maintenance of regions elevated in GC content) by biased gene conversion [57]. An alternative possibility is that base composition drives recombination rate variation and it is conservation of GC content, or GC-rich motifs [109], that results in the appearance of recombination rate conservation.

Further, the influence of Hill-Robertson interference on the accumulation of mildly deleterious substitutions is not expected to decrease linearly with an increase of the recombination rate [37, 110]. In this context, it is possible that the recombination rate is too high in most regions of the chicken and the zebra finch genome to lead to measurable variation in the efficiency of selection. This would somewhat contradict the observation of very low recombination rates in the sub-centromeric region of the larger zebra finch chromosomes [32]. However, these recombination deserts can have a high effective number of recombination events given a very large population size, as is observed for natural zebra finch populations [111]. In general, it may very well be that the effective population sizes of ancestral passerines have been higher than that of other (larger) birds.

Alignments

We downloaded protein-coding sequences from the chicken (G. gallus, WASHUC2), zebra finch (T. guttata, TaeGut3.2.4), green lizard (Anolis carolinensis, ANOCAR1), short-tailed opossum (Monodelphis domestica, MonDom5), platypus (Ornithorhynchus anatinus, OANA5), mouse (Mus musculus, NCBIM37) and human (Homo sapiens, NCBI36) genome assemblies through biomart [112] in Ensembl version 55. In order to identify 1:1 orthologs between zebra finch and each of the other species, we used a reciprocal Blast best hit approach as implemented in Inparanoid3.0 [113]. Codon-based pairwise alignments from the corresponding protein sequences were made using MUSCLE3.7 [114]. We used Gblocks 0.91b [115] to eliminate poorly aligned positions. In total, our analysis was based on 8,384 genes.

Estimates of substitution rates

Statistical models

We used bivariate and partial correlations to analyze the relationship between ω and recombination rate separately in chicken and zebra finch. The sex-average recombination rate for 1-Mb windows was obtained for chicken from Groenen et al. [31] and for zebra finch from Backström et al. [32]. Partial correlations controlled for GC content and the amount of coding sequence within each window individually and in combination. Similarly, we used bivariate and partial correlation (controlling for GC content) to study the association between divergence at fourfold degenerate sites from 1-Mb windows in different bird species. Since the windows were not identical between zebra finch and chicken, we estimated the correlations separately for zebra finch-chicken, and for chicken-zebra finch. The similarity in the results shows that the analysis is not susceptible to the exact location of windows. When controlling for GC content in correlations between zebra finch/chicken and the ancestral bird linage, we used the average GC content of both chicken and zebra finch as an estimate of GC content in the ancestral lineage.

Identification of candidate genes for adaptive evolution

Gene Ontology analysis

To test for overrepresentation of biological processes, molecular functions and cellular components among positively selected or rapidly evolving genes, we performed GO analysis using GoStat [121]. We downloaded GO annotations for chicken, human and mouse from Biomart. The analysis was based on Fisher's exact test between two lists of genes, that is, PS genes and a reference list of all analyzed genes. Multiple testing corrected significance values were based on Benjamini and Hochberg [122] correction (adjusted P < 0.1), included with the GoStat software.

Analysis of neurological genes

The OMIM database [123] was searched on 29 March 2009, using three different search phrases and a search limit set for 'prefix star' (that is, to find only OMIM terms associated with a known gene sequence). One search was on the term 'learning'. To search for genes related to neurogenesis, we used this phrase: [(stem cell AND neur*) OR neurogen*]. To search for genes related to neurodegeneration, we used 'neurogen*'. Human gene IDs in OMIM were cross-referenced and corrected or completed as needed against the HGNC database, and used to retrieve Ensembl gene IDs for human and zebra finch orthologs (Ensembl 53) via Biomart.