Birds are the most species-rich class of tetrapod vertebrates and have wide relevance across many research fields. We explored bird macroevolution using full genomes from 48 avian species representing all major extant clades. The avian genome is principally characterized by its constrained size, which predominantly arose because of lineage-specific erosion of repetitive elements, large segmental deletions, and gene loss. Avian genomes furthermore show a remarkably high degree of evolutionary stasis at the levels of nucleotide sequence, gene synteny, and chromosomal structure. Despite this pattern of conservation, we detected many non-neutral evolutionary changes in protein-coding genes and noncoding regions. These analyses reveal that pan-avian genomic diversity covaries with adaptations to different lifestyles and convergent evolution of traits.

With ~10,500 living species (1), birds are the most species-rich class of tetrapod vertebrates. Birds originated from a theropod lineage more than 150 million years ago during the Jurassic and are the only extant descendants of dinosaurs (2, 3). The earliest diversification of extant birds (Neornithes) occurred during the Cretaceous period. However, the Neoaves, the most diverse avian clade, later underwent a rapid global expansion and radiation after a mass extinction event ~66 million years ago near the Cretaceous-Paleogene (K-Pg) boundary (4, 5). As a result, the extant avian lineages exhibit extremely diverse morphologies and rates of diversification. Given the nearly complete global inventory of avian species, and the immense collected amount of distributional and biological data, birds are widely used as models for investigating evolutionary and ecological questions (6, 7). The chicken (Gallus gallus), zebra finch (Taeniopygia guttata), and pigeon (rock dove) (Columba livia) are also important model organisms in disciplines such as neuroscience and developmental biology (8). In addition, birds are widely used for global conservation priorities (9) and are culturally important to human societies. A number of avian species have been domesticated and are economically important. Farmed and wild water birds are key players in the global spread of pathogens, such as avian influenza virus (10).

Despite the need to better understand avian genomics, annotated avian genomic data was previously available for only a few species: the domestic chicken, domestic turkey (Meleagris gallopavo) and zebra finch (11–13), together with a few others only published recently (14–16). To build an understanding of the genetic complexity of birds and to investigate links between their genomic variation and phenotypic diversity, we collected and compared genome sequences of these and other avian species (48 species total), representing all 32 neognath and two of the five palaeognath orders (Fig. 1) (17), thus representing nearly all of the major clades of living birds (5).

Results

Sequencing, assembly, and annotation We used a whole-genome shotgun strategy to generate genome sequences of 45 new avian species (18), including two species representing two orders within the infraclass Paleognathae [common ostrich (Struthio camelus) and white-throated tinamou (Tinamus guttatus)], the other order within Galloanserae [Peking duck (Anas platyrhynchos)], and 41 species representing 30 neoavian orders (table S1) (19). In combination with the three previously published avian genomes (11–13), the genome assemblies cover 92% (34 of 37) of all avian orders (the three missing orders belong to the Paleognathae) (17). With the exception of the budgerigar (Melopsittacus undulatus), which was assembled through a multiplatform (Illumina/GS-FLX/PacBio) approach (20), all other new genomes were sequenced and assembled with Illumina (San Diego, CA) short reads (Fig. 1) (18). For 20 species, we produced high (>50×) coverage sequences from multiple libraries, with a gradient of insert sizes and built full-genome assemblies. For the remaining 25 species, we generated low (~30×) coverage data from two insert-size libraries and built less complete but still sufficient assemblies for comparative genome analyses. These de novo (18) genome assemblies ranged from 1.05 to 1.26 Gb, which is consistent with estimated cytology-based genome sizes (21), suggesting near complete genome coverage for all species. Scaffold N50 sizes for high-coverage genomes ranged from 1.2 to 6.9 Mb, whereas those for lower-coverage genomes were ~48 kb on average (table S2). The genomes of the ostrich and budgerigar were further assembled with optical maps, increasing their scaffold N50 sizes to 17.7 and 13.8 Mb, respectively (20, 22). We annotated the protein-coding sequences using a homology-based method for all genomes, aided by transcriptome sequencing for some species (18). To avoid systematic biases related to the use of different methods in annotations of previously published avian genomes, we created a uniform reference gene set that included all genes from the chicken, zebra finch, and human (23). This database was used to predict protein gene models in all avian genomes and American alligator (Alligator mississippiensis) (24). All high-coverage genomes were predicted to contain ~15,000 to 16,000 transposable element-free protein-coding genes [table S3 and annotation files in (19)], similar to the chicken genome (~15,000). Despite the fragmented nature of the low-coverage genomes leading to ~3000 genes likely missing or partially annotated, it was still possible to predict 70 to 80% of the entire catalog of avian genes.

Broad patterns of avian genome evolution Although many fishes and some amphibians have smaller genomes than birds, among amniotes, birds have the smallest (21). The genomes of mammals and nonavian reptiles typically range from 1.0 to 8.2 Gb, whereas avian genomes range from 0.91 in the black-chinned hummingbird (Archilochus alexanderi) to a little over 1.3 Gb in the common ostrich (21). A number of hypotheses have been proposed for the smaller avian genome size (25–28). Here, we document key events that have likely contributed to this smaller genome size. The proliferation and loss of transposable elements (TEs) may drive vertebrate genome size evolution (29–31). Consistent with the zebra finch and galliformes genomes (11–13, 32), almost all avian genomes contained lower levels of repeat elements (~4 to 10% of each genome) (table S4) than in other tetrapod vertebrates (for example, 34 to 52% in mammals) (33). The sole outlier was the downy woodpecker (Picoides pubescens), with TEs representing ~22% of the genome, derived mainly from species-specific expansion of LINE (long interspersed elements) type CR1 (chicken repeat 1) transposons (fig. S1). In contrast, the average total length of SINEs (short interspersed elements) in birds has been reduced to ~1.3 Mb, which is ~10 to 27 times less than in other reptiles [12.6 Mb in alligator; 34.9 Mb in green sea turtle (Chelonia mydas)], suggesting that a deficiency of SINEs occurred in the common ancestor of birds. We compared the average size of genomic elements of birds with 24 mammalian and the three nonavian reptile genomes. Avian protein-coding genes were on average 50 and 27% shorter than the mammalian and reptilian genes, respectively (Fig. 2A). This reduction is largely due to the shortening of introns and reduced intergenic distances that resulted in an increased gene density (Fig. 2A). Such genomic contraction has also evolved convergently in bats (fig. S11), the only flying mammalian group. The condensed genomes may represent an adaptation tied to rapid gene regulation required during powered flight (34, 35). Fig. 2 Genome reduction and conservation in birds. (A) Comparison of average size of introns, exons, and intergenic regions within avian, reptilian, and mammalian genomes. (B) Synteny plot and large segmental deletions between green anole chromosome 2 and multiple chicken chromosomes. Colored bars and lines indicate homologous blocks between two species; black bars indicate location of large avian-specific segmental deletions, which are enriched at the breakpoints of interchromosome rearrangements. (Bottom) An example of a large segmental deletion in birds (represented by ostrich genes). Homologous genes annotated in each species are shown in small boxes. The color spectrum represents the percent identity of homologous genes with the green anole. (C) Distribution of gene synteny percentages identified for phylogenetically independent species pairs of various divergence ages. Dots indicate the percentage of genes remaining in a syntenic block in pairwise comparisons between two avian or mammalian species. Box plots indicate that the overall distributions of the synteny percentages in birds and mammals are different (P value was calculated by using Wilcoxon rank sum test with phylogenetically independent species pairs). (D) Chromosomal organization of the α- and β-globin gene clusters in representative avian and mammalian taxa. These genes encode the α- and β-type subunits of tetrameric (α2β2) hemoglobin isoforms that are expressed at different ontogenetic stages. In the case of the α-like globin genes, birds and mammals share orthologous copies of the αD- and αA-globin genes. Likewise, the avian π-globin and the mammalian ζ-globin genes are 1:1 orthologs. In contrast, the genes in the avian and mammalian β-globin gene clusters are derived from independent duplications of one or more β-like globin genes that were inherited from the common ancestor of tetrapod vertebrates (90, 91). To further investigate whether avian genome size reduction is due to a lineage-specific reduction in the common avian ancestor of birds or expansion in other vertebrates (36), we performed ancestral state reconstructions of small [<100 base pairs (bp)] deletion events across an alignment of four representative well-assembled avian and three reptile genomes (18) and found that the avian ancestral lineage experienced the largest number of small deletion events—about twice the number in the common ancestor of birds and crocodiles (fig. S12). In contrast, many fewer small deletion events occurred in modern avian lineages (fig. S12). We next created a gene synteny map between the highest-quality assembled avian genome (ostrich) and other reptile genomes to document lineage-specific events of large segmental deletions (18). We detected 118 syntenic blocks, spanning a total of 58 Mb, that are present in alligator and turtle genomes but lost in all birds (table S8). In contrast, ~8x and ~5x fewer syntentic blocks were missing in alligator (14 blocks, 9 Mb) and turtle (27 blocks, 8 Mb) relative to green anole, respectively, confirming the polarity of genome size reduction in birds (table S8). The large segmental losses in birds were skewed to losses from chr2 and chr6 of the green anole (fig. S13). Two of the green anole’s 12 pairs of microchromosomes, LGd and LGf, were completely missing in birds, with no homologous genes found within the avian genomes. Most of these lost segments were located at the ends of chromosomes or close to the centrosomes (fig. S13). Furthermore, lost segments were enriched at apparent breakpoints of the avian microchromosomes (Fig. 2B and fig. S13). These findings imply that the large segmental losses may be a consequence of chromosomal fragmentation events in the common ancestor of birds giving rise to additional microchromosomes in modern birds. The large segmental deletions in birds contain at least 1241 functional protein-coding genes (table S9), with each lost segment containing at least five contiguous genes. The largest region lost in birds was a 2.1-Mb segment of the green anole chr2, which contains 28 protein-coding genes (Fig. 2B). Overall, at least 7% of the green anole macrochromosomal genes were lost through segmental deletions in birds. Although gene loss is a common evolutionary process, this massive level of segmental deletion has not been previously observed in vertebrates. Over 77% of the 1241 genes present in the large segmentally deleted regions have at least one additional paralog in the green anole genome, a level higher than the overall percentage of genes with paralogs in the green anole genome or avian genomes (both at ~70%). This suggests that birds may have undergone functional compensation in their paralogous gene copies, reducing selection against the loss of these segmental regions. We predict that the loss of functions associated with many genes in the avian ancestor may have had a profound influence on avian-specific traits (table S11).

Conservative mode of genome evolution With ~2/3 of avian species possessing ~30 pairs of microchromosomes, the avian karyotype appears to be distinctly conserved because this phenotype is not a general feature of any other vertebrate group studied to date (37). We assessed the rates of avian chromosomal evolution among the 21 more fully assembled genomes (scaffold N50 > 1 Mb) (table S2) (18). From the alignment of chicken with the other 20 avian genomes, plus green anole and Boa constrictor (38), we identified homologous synteny blocks (HSBs) and 1746 evolutionary breakpoint regions (EBRs) in different avian lineages and then estimated the expected number of EBRs (18) and the rates of genomic rearrangements, using a phylogenetic total evidence nucleotide tree (TENT) as a guide (5). We excluded the turkey genome after detecting an unusually high fraction of small lineage-specific rearrangements, suggesting a high number of local misassemblies. Of the 18 remaining non–Sanger-sequenced genomes (table S2), the estimated rate of chimeric scaffolds that could lead to false EBRs was ~6% (39). The average rate of rearrangements in birds is ~1.25 EBRs per million years; however, bursts of genomic reorganization occurred in several avian lineages (fig. S15). For example, the origin of Neognathae was accompanied by an elevated rate of chromosome rearrangements (~2.87 EBRs per million years). Intriguingly, all vocal learning species [zebra finch, medium-ground finch (Geospiza fortis), American crow (Corvus brachyrhynchos), budgerigar, and Anna’s hummingbird (Calypte anna)] had significantly higher rates of rearrangements than those of close vocal nonlearning relatives [golden-collared manakin (Manacus vitellinus), peregrine falcon (Falco peregrinus) and chimney swift (Chaetura pelagica)] [phylogenetic analysis of variance F statistic (F) = 5.78, P = 0.0499] and even higher relative to all vocal nonlearning species (F = 15.03, P = 0.004). This may be related to the larger radiations these clades experienced relative to most other bird groups. However, the golden-collared manakin, which belongs to suboscines (vocal nonlearners) that have undergone a larger radiation than parrots and hummingbirds, has a low rearrangement rate. We next compared microsynteny (local gene arrangements), which is more robust and accurate than macrosynteny analyses for draft assemblies (18). We compared with eutherian mammals, which are approximately the same evolutionary age as Neoaves and whose genome assemblies are of similar quality. We examined the fraction of orthologous genes identified from each pair of two-avian/mammalian genomes, on the basis of syntenic and best reciprocal blast matches (18). Birds have a significantly higher percentage of synteny-defined orthologous genes than that of mammals (Fig. 2C). The fraction of genes retained in syntenic blocks in any pairwise comparison was linearly related with evolutionary time, by which the overall level of genome shuffling in birds was lower than in mammals over the past ~100 million years (Fig. 2C). This suggests a higher level of constraint on maintaining gene synteny in birds relative to mammals. The apparent stasis in avian chromosome evolution suggests that birds may have experienced relatively low rates of gene gain and loss in multigene families. We examined the intensively studied gene families that encode the various α- and β-type subunits of hemoglobin, the tetrameric protein responsible for blood oxygen transport in jawed vertebrates (40). In amniotes, the α- and β-globin gene families are located on different chromosomes (40) and experienced high rates of gene turnover because of lineage-specific duplication and deletion events (41). In birds, the size and membership composition of the globin gene families have remained remarkably constant during ~100 million years of evolution, with most examined species retaining an identical complement (Fig. 2D). Estimated gene turnover rates (λ) of α- and β-globin gene families were over twofold higher in mammals than birds (λ = 0.0023 versus 0.0011, respectively). Much of the variation in the avian α-globin gene family was attributable to multiple independent inactivations of the αD-globin gene (Fig. 2D), which encodes the α-chain subunit of a hemoglobin isoform (HbD) expressed in both embryonic and definitive erythrocytes (42). Because of uniform and consistent differences in oxygen-binding properties between HbD and the major adult-expressed hemoglobin isoform, HbA (which incorporates products of αA-globin) (42), the inactivations of αD-globin likely contribute to variation in blood-oxygen affinity, which has important consequences for circulatory oxygen transport and aerobic energy metabolism. Overall, the globin gene families illustrate a general pattern of evolutionary stasis in birds relative to mammals. Genomic nucleotide substitution rates vary across species and are determined through both neutral and adaptive evolutionary processes (43). We found that the overall pan-genomic background substitution rate in birds (~1.9 × 10–3 substitutions per site per million years) was lower than in mammals (~2.7 × 10–3 substitution per site per million years) (Fig. 3A). However, the substitution rate estimates also exhibited interordinal variation among birds (Fig. 3A). There was a positive correlation between the substitution rate and the number of species per order [coefficient of determination (R2) = 0.21, P = 0.01, Pearson’s test with phylogenetically independent contrasts] (Fig. 3B and fig. S19), evidencing an association with rates of macroevolution (44). For example, Passeriformes, the most diverse avian order, exhibited the highest evolutionary rate (~3.3 × 10–3 substitutions per site per million years), almost two times the average of Neoaves (~2 × 10–3 substitutions per site per million years, Fig. 3A). Landbirds exhibited an average higher substitution rate than that of waterbirds (landbirds, ~2.2 × 10–3 substitutions per site per million years; waterbirds, ~1.6 × 10–3 substitutions per site per million years), which is consistent with the observation that landbirds have greater net diversification rates than those of waterbirds (7). Among the landbirds, the predatory lineages exhibited slower rates of evolution (~1.6 × 10–3 substitutions per site per million years), similar to that of waterbirds. Moreover, the three vocal learning landbird lineages (parrots, songbirds, and hummingbirds) are evolving faster than are nonvocal learners (Fig. 3A). Overall, our analyses indicate that genome-wide variation in rates of substitution is a consequence of the avian radiation into a wide range of niches and associated phenotypic changes. Fig. 3 Evolutionary rate and selection constraints. (A) Substitution rate in each lineage was estimated by the comparison of fourfold degenerate (4d) sites in coding regions, in units of substitutions per site per million years. Waterbirds and landbirds are defined in (5). (B) Correlation between average substitution rates and number of species within different avian orders. Divergence times were estimates from (5). The fit line was derived from least square regression analysis, and the confidence interval was estimated by “stat_smooth” in R. The units of the x axis are numbers of substitutions per site per million years. The correlation figure with phylogenetically independent contrasts is provided in the supplementary materials. (C) Density map for comparison of conservation levels between pan-avian and pan-mammalian genomes, on the basis of the homologous genomic regions between birds and mammals. Conservation levels were quantified by means of PhastCons basewise conservation scores. (D) HCEs found in both mammalian and avian genomes (smaller pie piece) and those that are avian-specific (larger pie piece). (E) MID1 contains abundant avian-specific HCEs in the upstream and downstream regulatory regions. Many regulatory motif elements are identified in these avian-specific HCEs. Cons., conservation level.

Selective constraints on functional elements Conservation of DNA sequences across distantly related species reflects functional constraints (45). A direct comparison of 100-Mb orthologous genomic regions revealed more regions evolving slower than the neutral rate among birds (Fig. 3C) than mammals (46), which is consistent with the slower rate of avian mitochondrial sequence evolution (47). We predicted 3.2 million highly conserved elements (HCEs) at a resolution of 10 bp or greater spanning on average 7.5% of the avian genome, suggesting a strong functional constraint in avian genomes. Functional annotations revealed that ~12.6% of these HCEs were associated with protein-coding genes, whereas the majority of the remaining HCEs were located in intron and intergenic regions (Fig. 3, D and E). These HCEs enabled us to identify 717 new protein-coding exons and 137 new protein-coding genes, with 77% of the latter supported by the deep transcriptome data (table S17). Deep transcriptome sequencing also enabled us to annotate 5879 candidate long noncoding RNA (lncRNA) genes, of which 220 overlapped HCEs with a coverage ratio of >50% (table S18) (18). Because HCEs may have different functions in different lineages, we separated the HCEs into two categories: bird-specific and amniote HCEs (shared by birds and mammals). Among the bird-specific HCEs, we identified 13 protein-coding genes that were highly conserved in birds but divergent in mammals (table S19). One of the most conserved was the sperm adhesion gene, SPAM1, which mediates sperm binding to the egg coat (48). This gene, however, was under positive selection driven by sperm competition in mammalian species (49). Noncoding HCEs play important roles in the regulation of gene expression (50); thus, we compared the transcription factor binding sites in the ENCODE project (51) with the HCEs and found that the avian-specific HCEs are significantly associated with transcription factors functioning in metabolism (table S20), whereas amniote core HCEs are enriched with transcription factors functioning in signal regulation, stimulus responses, and development (table S21). To investigate evolutionary constraints on gene regions, we calculated dN/dS [the ratio of the number of nonsynonymous substitutions per nonsynonymous site (dN) to the number of synonymous substitutions per synonymous site (dS)] for 8295 high-quality orthologs. Consistent with the fast-Z sex chromosome hypothesis (52), the evolutionary rate of Z-linked genes was significantly higher than autosome genes (Fig. 4A). This is most likely driven by the reduction of effective population size (N e ) of Z-linked genes—because the N e of Z chromosome is only 3/4 of that of autosomes—as well as by male sexual selection (52). Furthermore, consistent with the fast-macro hypothesis, the overall rate of macrochromosomal genic evolution is higher than that of microchromosomes (Fig. 4A), which is probably due to differences in the recombination rates and genic densities between macro- and microchromosomes in birds (53). Fig. 4 Selection constraints on genes. (A) Box plot for the distribution of dN/dS values of genes on avian macrochromosomes, microchromosomes, and the Z chromosome. P values were calculated with Wilcoxon rank sum tests. (B) GO categories in Neoaves, Galloanserae, and Palaeognathae showing clade-specific rapid evolutionary rates. Red bars, P value of significance; blue bars, number of genes in each GO. We also examined the dN/dS ratio of each avian Gene Ontology (GO) category for comparison with mammals and within birds. Those involved in development (such as spinal cord development and bone resorption) are evolving faster in birds, and those involved in the brain function (such as synapse assembly, synaptic vesicle transport, and neural crest cell migration) are evolving faster in mammals (tables S23 and S24). Genes involved in oxidoreductase activity were relatively rapidly evolving in the Palaeognathae clade that contains the flightless ratites (Fig. 4B and table S25). The fast evolving GOs in the Galloanserae participate in regulatory functions (Fig. 4B and table S26). In Neoaves, genes involved in microtubule-based processes were the fastest evolving (Fig. 4B and table S27). We speculate that these differences could be caused by relaxed selective constraints or positive selection in different lineages.

Genotype-phenotype convergent associations: Evolution of vocal learning With the availability of genomes representing all major modern avian lineages and their revised phylogenetic relationships (5), it becomes possible to conduct genome-wide association studies across species with convergent traits. We focused on vocal learning, which given our phylogenetic analyses is inferred as having evolved independently, either twice, in hummingbirds and the common ancestor of songbirds and parrots, or three times (5, 54). All three groups have specialized song-learning forebrain circuits (song nuclei) not found in vocal nonlearners (Fig. 5A) (55). Fig. 5 Convergent molecular changes among vocal learning birds. (A) Songbird brain diagram showing the specialized forebrain song-learning nuclei (yellow) that controls the production (HVC and RA) and acquisition (LMAN and Area X) of learned song (55). Gray arrows indicate connections between brain regions; red and blue (thick) arrows indicate relative numbers of genes with increased or decreased specialized expression in zebra finch song nuclei and with convergent accelerated coding sequences (left numbers of 66 total) or convergent amino acid substitutions (right numbers of 6). Genes expressed in more than one song nucleus are counted multiple times. RA, robust nucleus of the arcopallium; LMAN, lateral magnocellular nucleus of the anterior nidopallium; Area X, Area X of the striatum; and HVC, a letter-based name. (B) Classification of vocal learner-specific accelerated elements, compared with the background alignment of 15 avian species. Analyses of 7909 orthologous protein-coding genes with available amino acid sites in all three vocal-learning and control vocal nonlearning groups revealed convergent accelerated dN/dS for 227 genes in vocal learners (table S28). Of these, 73% (165) were expressed in the songbird brain (physically cloned mRNAs), and of these, 92% (151) were expressed in adult song-learning nuclei, which is much higher than the expected 60% of brain genes expressed in song nuclei (56). About 20% (33) were regulated by singing, which is twice the expected 10% (56). In addition, 41% of the song nuclei accelerated genes showed differential expression among song nuclei [expected 20% (56)], and 0.7 to 9% [0.7 to 4.3% expected (57)] showed specialized expression compared with the surrounding brain regions (table S28) (58). GO analyses of the accelerated differentially expressed song nuclei genes revealed 30 significant functionally enriched gene sets, which clustered into four major categories, including neural connectivity, brain development, and neural metabolism (fig. S25). For an independent measure of convergence, we developed an approach that scans for single amino acid substitutions common to species with a shared trait, controlling for phylogenetic relationships (18). Of the 7909 genes, 38 had one to two amino acid substitutions present only in vocal learners (table S31). At least 66% of these were expressed in the songbird brain, including in the song nuclei [58%; 20% expected (56)]. Two genes (GDPD4 and KIAA1919) showed convergent accelerated evolution on the amino acid sites specific to vocal learners (table S31). To identify accelerated evolution in noncoding sequences in vocal learners, we scanned the genome alignment using phyloP (18, 59). We used a more limited sampling of vocal nonlearning species closely related to the vocal learners (table S32) because of the relatively faster evolutionary rate of noncoding regions. We scanned the entire genome alignment and found 822 accelerated genomic elements specifically shared by all three vocal learning groups (table S33). These convergent elements were skewed to intergenic regions in vocal learners relative to the background average accelerated elements across species (Fisher’s exact test, P < 2.2 × 10–16) (Fig. 5B). Of these elements, 332 were associated with 278 genes (within 10 kb 5′ or 3′ of the nearest gene), of which a high proportion (76%) was expressed in the brain; almost all of those (94%, 198 genes) expressed in one or more song nuclei, 20% were regulated by singing (10% expected), 51% (20% expected) showed differential expression among song nuclei, and 2 to 15% [0.7 to 4.3% expected, based on (56)] had specialized expression relative to the surrounding brain regions, including the FoxP1 gene involved in speech (table S34 and figs. S27 to S32). Overall, these analyses show a 2- to 3.5-fold enrichment of accelerated evolution in regulatory regions of genes differentially expressed in vocal learning brain regions. In contrast, there was very little overlap (2.5%) of genes with convergent accelerated noncoding changes and convergent accelerated amino acid changes, indicating two independent targets of selection for convergent evolution.