Dense sampling of bird diversity increases power of comparative genomics

Abstract

Whole-genome sequencing projects are increasingly populating the tree of life and characterizing biodiversity^1-4. Sparse taxon sampling has previously been proposed to confound phylogenetic inference⁵, and captures only a fraction of the genomic diversity. Here we report a substantial step towards the dense representation of avian phylogenetic and molecular diversity, by analysing 363 genomes from 92.4% of bird families-including 267 newly sequenced genomes produced for phase II of the Bird 10,000 Genomes (B10K) Project. We use this comparative genome dataset in combination with a pipeline that leverages a reference-free whole-genome alignment to identify orthologous regions in greater numbers than has previously been possible and to recognize genomic novelties in particular bird lineages. The densely sampled alignment provides a single-base-pair map of selection, has more than doubled the fraction of bases that are confidently predicted to be under conservation and reveals extensive patterns of weak selection in predominantly non-coding DNA. Our results demonstrate that increasing the diversity of genomes used in comparative studies can reveal more shared and lineage-specific variation, and improve the investigation of genomic characteristics. We anticipate that this genomic resource will offer new perspectives on evolutionary processes in cross-species comparative analyses and assist in efforts to conserve species.

Fig. 1. Newly sequenced genomes densely cover the bird tree of life.

The 10,135 bird species^, are shown on a draft phylogeny that synthesizes taxonomic and phylogenetic information (Supplementary Data). In total, 363 species, covering 92.4% of all families, now have at least 1 genome assembly per sequenced family (purple branches). The grey arc marks the diverse Passeriformes radiation, with 6,063 species, of which 173 species have genome assemblies now. Chicken (*) and zebra finch (**) are marked for orientation. Paintings illustrate examples of sequenced species.

Fig. 2. Improved orthologue distinction and detection of lineage-specific sequences.

a, Incorporating synteny in the orthologue assignment pipeline resolves complex cases of orthology. The growth hormone gene (GH) has one copy in chicken and two copies in Passeriformes (exemplified by zebra finch and Atlantic canary). On the basis of the conserved synteny of the GH_L in Passeriformes with GH of chicken, the pipeline recognized GH_L as the ancestral copy—despite high similarity to the other copy. b, The whole-genome alignment allows detecting lineage-specific sequences. For orders with more than one sequenced representative, lineage-specific sequences are those present in the reconstructed ancestral genome but absent in other lineages. Colours denote higher-level taxonomic groupings. c, A novel gene in Passeriformes. Phylogeny based on the B10K Project phase I plotted with synteny of a putative lineage-specific gene (DNAJC15L) and its surrounding genes. DNAJC15L is found in 131 out of 173 sequenced Passeriformes and their reconstructed ancestral genome, but is not found in non-Passeriformes. MRCA, most-recent common ancestor.

Fig. 3. Denser phylogenomic sequencing increases the power to detect selective constraints.

Results are shown from 3 alignments for 53 birds, 77 vertebrates, and 363 birds. a, Proportion of alignment columns labelled as conserved. The cumulative portion of the genome with a conserved call is shown, starting from the column with the smallest P value and proceeding to the columns with the highest P values. The dotted lines show the path after hitting the false-discovery rate (FDR) P value cutoff of 0.05, below which calls are significant (marked by arrows). b, Histograms of the rate of alignment columns evolving slower relative to the neutral rate (labelled 1.00). Coloured areas indicate significantly conserved columns, and light grey areas indicate non-significantly conserved columns. A rate of zero contains a relatively high proportion of recent insertions present in only a few species; there is limited statistical power to classify such insertions. c, Proportion of various functional regions of the chicken genome that contain single-bp conserved elements in the large alignment compared to alignments with fewer species. UTR, untranslated region. d, An example of a MAF::NFE2 motif overlaid on one of its predicted binding sites demonstrates the high resolution of our conserved site predictions and the increased power to predict conservation in the larger alignment.

Extended Data Fig. 1. Sampling and processing of the 363 genomes.

a, Sources of the 363 genomes. Each genome is a square; colour indicates the data source. Newly published genomes from the B10K Project phase II are red; unpublished genomes contributed by external labs are yellow; published genomes from phase I are orange; genomes contributed by the community that have since been published are dark blue; and other genomes available on NCBI are light blue. b, Map of geographical origin of the 281 bird samples for which geographical coordinates are available. c, Summary of the species confirmation of 236 B10K Project newly sequenced species. The downward arrows are excluded genomes. d, Summary of mitochondrial genome assembly and annotation for 336 species. The downward arrows are excluded mitochondrial genomes.

Extended Data Fig. 2. Distribution of transposable elements.

a, Percentage of the genome that is a transposable element (TE). Box plots are shown for groups with at least three sequenced species. b, Per cent base pairs of the genome that are long interspersed nuclear elements (LINEs), grouped by orders. Box plots are shown for groups with at least three sequenced species. c, S.d. of the transposable element content for orders with at least three sequenced species. d, S.d. of the per cent LINE content for orders with at least three sequenced species. e, Ancestral state reconstruction of total transposable elements. The branch colour from blue to red indicates an increase in transposable elements. Two orders with noticeable patterns—Piciformes and Bucerotiformes—are labelled on the tree. A zoomable figure with labels for all terminals is available at www.doi.org/10.17632/fnpwzj37gw.

Extended Data Fig. 3. Patterns of the presence and absence of 5 visual opsins in 363 bird species.

This figure shows patterns for the visual opsins encoded by RH1, RH2, OPN1sw1, OPN1sw2 and OPN1lw. Colours correspond to five annotated states of opsin sequences. A zoomable figure with labels for all terminals is available at www.doi.org/10.17632/fnpwzj37gw.

Extended Data Fig. 4. GC content and codon use.

a, Principal component analysis (PCA) of GC content in the coding regions of orthologues with conserved synteny with chicken for 340 bird species, including 164 Passeriformes species. b, Correspondence analysis of RSCU for all 363 birds. The primary and secondary axes account for 78.18% and 14.82% of the total variation, respectively. c, The distribution of codons on the same two axes as shown in b, with each codon coloured according to its ending nucleotide. This showed that the axis-1 score of a species is primarily determined by differences in frequencies of codons ending in G, C, A or T. d, RSCU analysis of 59 codons across avian genomes (n = 363 biologically independent species for each box plot). The horizontal lines indicate thresholds of under-represented codons (<0.6, blue box plots), average representation (1.0, white box plots) and over-represented codons (>1.6, orange box plots). e, Pearson correlation between GC content of the third codon position and the primary axis in b, colour-coded to distinguish Passeriformes and non-Passeriformes. The strong correlation (R² = 0.9, P = 4.1 × 10⁻¹⁸⁴) indicates that the frequencies of codons ending in G or C is the main driver of the codon bias in Passeriformes. f, Comparison of the mean Nc values between the Passeriformes and other species for orthologues with conserved synteny with chicken (Supplementary Table 12). Each dot represents the mean Nc value of an orthologue in the Passeriformes and other species, respectively. Orthologues with at least 20 individuals in both the Passeriformes and the non-Passeriformes were included in this analysis.

Extended Data Fig. 5. Overview of the pipelines for identifying genomic regions.

a, Assignment of orthologous protein-coding regions. All pairwise relationships between homologous regions obtained from the Cactus alignment (4 species shown here in different colours) were used to construct the homologous groups across all 363 birds. Using chicken as the reference, we further generated a table containing homologues with conserved synteny to chicken. b, Annotation of conserved orthologous intron regions on the basis of Cactus whole-genome alignments. The credible intron fragments in chicken were picked out after filtering out regions mapped by RNA sequences, and chicken-specific or repetitive regions. Orthologous relationships of intron fragments were detected on the basis of the aligned Cactus hits and the orthologues with conserved synteny with chicken. The non-intron regions of each bird in the alignments were masked as gaps.

Extended Data Fig. 6. Gene tree for copies of the growth hormone gene GH.

The tree was generated by maximum likelihood phylogenetic analysis of avian GH gene copies. Only nodes with >80 bootstrap are annotated as dots; the larger the dot, the higher the bootstrap. All Passeriformes sequences are clustered in a single clade and there are two sister gene clades within Passeriformes, corresponding to the GH_S gene copy (blue) and the GH_L gene copy (orange). Twelve species with only one copy are indicated by green stars. A zoomable figure with labels for all terminals and the tree file is available at www.doi.org/10.17632/fnpwzj37gw.

Extended Data Fig. 7. Identification of lineage-specific sequences.

a, An example of a 36-bp insertion (red) identified by Cactus in the southern cassowary (Casuarius casuarius) compared to the Okarito brown kiwi (Apteryx rowi) (both in Palaeognathae) with mapped sequence reads shown as lines. b, Proportion of lineage-specific sequence for each order correlated with the distance from parent node to MRCA node (branch length). c, Presence and absence of the DNAJC15-like gene (DNAJC15L), and its surrounding genes, in all 363 birds. Upstream: KLHL1 and DACH1; downstream: MZT1, BORA, RRP44, PIBF1 and KLF5. The state is shown for each bird in three ways: multiple copies (filled shapes), one copy (empty shapes) and no gene (blank). Passeriformes are highlighted in red. A zoomable figure with labels for all terminals is available at www.doi.org/10.17632/fnpwzj37gw. d, Exon fusion patterns of the DNAJC15-like gene (DNAJC15L) in three Passeriformes, compared to exon structure of the ancestral DNAJC15. For L. aspasia, gene models for the ancestral and novel copy are shown. The structure of the ancestral copy is highly conserved across all bird species with five introns. The Passeriformes-specific copy has no intron or newly derived minor intron and includes a poly-(A) at the 5′ end, which implies that this new gene was derived from retroduplication of DNAJC15.

Extended Data Fig. 8. The evolution of songbirds was associated with the loss of the cornulin gene.

a, Presence and absence of the cornulin gene (CRNN) and its surrounding genes (EDDM and S100A11) in all 363 birds. Branches are coloured as oscine Passeriformes (blue), non-oscine Passeriformes (green) and non-Passeriformes (black). The states of genes are shown in three ways: functional gene (filled box), pseudogene (empty box) and gene not found (blank). Genes were identified by Exonerate using phylogenetically diverse EDDM, CRNN and S100A11 sequences as queries. A zoomable figure with labels for all terminals is available at www.doi.org/10.17632/fnpwzj37gw. b, Hypothesis on the evolutionary loss of cornulin and the appearance of a fine-tuned extensibility of the oesophagus as a vocal tract filter in songbirds.

Extended Data Fig. 9. Acceleration and conservation scores.

Results are shown from 3 alignments for 53 birds, 77 vertebrates, and 363 birds. a, Acceleration (left) and conservation (right) within alignment columns on chicken. This panel is similar to Fig. 3a, but includes accelerated columns. b, Proportion of chicken functional regions covered by significantly accelerated or conserved sites. This panel is similar to Fig. 3c, but includes accelerated columns.

Extended Data Fig. 10. Distribution of acceleration and conservation scores.

a, Distribution of conservation and acceleration scores within different functional region types across alignments. Lines mark quartiles of the density estimates. b, Larger histogram of chicken column rates. This panel is similar to Fig. 3b, but includes accelerated columns ending at a rate of 10× the neutral rate. c, Difference in PhyloP scores (compared to original scores) after realignment with MAFFT for a random sample of significantly conserved sites. d, Comparison of the distribution of PhyloP scores across alignments. Scores indicate log-scaled probabilities of conservation (positive values) or acceleration (negative values) for each base in the genome. a and d show results from three alignments for 53 birds, 77 vertebrates and 363 birds.