The Genome of the "Great Speciator" Provides Insights into Bird Diversification.

Among birds, white-eyes (genus Zosterops) have diversified so extensively that Jared Diamond and Ernst Mayr referred to them as the "great speciator." The Zosterops lineage exhibits some of the fastest rates of species diversification among vertebrates, and its members are the most prolific passerine island colonizers. We present a high-quality genome assembly for the silvereye (Zosterops lateralis), a white-eye species consisting of several subspecies distributed across multiple islands. We investigate the genetic basis of rapid diversification in white-eyes by conducting genomic analyses at varying taxonomic levels. First, we compare the silvereye genome with those of birds from different families and searched for genomic features that may be unique to Zosterops. Second, we compare the genomes of different species of white-eyes from Lifou island (South Pacific), using whole genome resequencing and restriction site associated DNA. Third, we contrast the genomes of two subspecies of silvereye that differ in plumage color. In accordance with theory, we show that white-eyes have high rates of substitutions, gene duplication, and positive selection relative to other birds. Below genus level, we find that genomic differentiation accumulates rapidly and reveals contrasting demographic histories between sympatric species on Lifou, indicative of past interspecific interactions. Finally, we highlight genes possibly involved in color polymorphism between the subspecies of silvereye. By providing the first whole-genome sequence resources for white-eyes and by conducting analyses at different taxonomic levels, we provide genomic evidence underpinning this extraordinary bird radiation.

Introduction

One of the most striking biodiversity patterns is that certain groups of organisms have more species than others. However, the underlying causes for the uneven distribution of diversity are just being unraveled for a few taxa (Hughes and Eastwood 2006; Brawand et al. 2014). Among terrestrial vertebrates, birds are the most diverse with about 10,000 species (Del Hoyo et al. 2011), and within birds, diversity is also highly skewed, from just five species of kiwis (order Apterygiformes) to half of the entire bird species richness found within passerines (order Passeriformes). One of the most diverse passerine families is Zosteropidae, with over 120 species from 12 genera, of which the white-eyes (genus Zosterops) have the highest number of representatives. White-eyes are among the largest avian genera, with about 80 described species found throughout Africa, Australasia, and the Southern Pacific ocean (Van Balen 2008; Clements et al. 2014). Indeed, Zosterops is one of the fastest diversifying lineages of terrestrial vertebrates (Moyle et al. 2009; Valente et al. 2010; Jetz et al. 2012), hence the group has been considered one of the “great speciators” (Diamond et al. 1976).

White-eyes are exceptional island colonizers; they are the group of passerines that has colonized more islands worldwide, and have thus been a model of choice for studying rapid speciation, phenotypic evolution, and genetic changes following the colonization of new environments (Black 2010; Clegg 2010; Milá et al. 2010; Melo et al. 2011; Cox et al. 2014). In particular, one of the key outstanding questions regarding white-eye evolution is what is known as the paradox of the great speciator (Diamond et al. 1976): Despite the extraordinary dispersal ability of white-eyes, a substantial degree of differentiation generally occurs across islands, with a large proportion of species or subspecies being single island endemics (Mees 1957; Melo et al. 2011). For all these reasons, white-eyes are an ideal group of birds to disentangle the genomic underpinning of species diversification; however, genome-scale data for these taxa have been lacking (Bourgeois et al. 2013).

A particularly interesting species of white-eye to study microevolution and diversification is the silvereye, Zosterops lateralis. One of the most widely distributed species of Zosterops, the silvereye occurs in mainland Australia and has naturally colonized several islands and archipelagos in the South Pacific, including New Caledonia, Vanuatu, Fiji, Tasmania, and New Zealand. The silvereye displays extensive genetically based phenotypic variation across numerous islands (Clegg, Degnan, Kikkawa, et al. 2002; Clegg, Degnan, Moritz, et al. 2002) and is considered to be a species complex consisting of multiple subspecies and allopatric forms that vary in body size and other morphological characters (Clegg, Degnan, Moritz, et al. 2002). Although molecular phylogenetic studies have shed light on the patterns of microevolution in the silvereye (Clegg and Phillimore 2010), genome-wide resources have not been available.

Here, we present a high-quality reference genome of the silvereye and use it as a model to address hypotheses as to why white-eyes have diversified so dramatically. In combination with whole genome resequencing and restriction site associated DNA (RAD; Baird et al. 2008) data for other white-eye species, we examine the genetic underpinning of rapid diversification in the genus. We conducted genomic analyses at different taxonomic levels.

To investigate what differentiates the genome of Zosteropidae from those of other bird families, we use a comparative genomics approach. We can formulate several hypotheses as to why white-eyes have diversified so dramatically compared with other bird clades. For example, there have been several studies linking substitution rates, phenotypic change, and species diversity (Lanfear et al. 2010; Hugall and Stuart-Fox 2012); we would therefore expect white-eyes to have particularly high rates of DNA evolution. Genomic innovations and evolutionary flexibility may arise from gene gains or losses (Lynch 2000; Brawand et al. 2014). We may also find genes contributing to reproductive isolation or affecting dispersal abilities to evolve more frequently under positive selection in white-eyes relative to other avian taxa (Melo et al. 2011).

In order to study interspecific patterns of genomic differentiation within Zosterops, we compared the genomes of different species of white-eyes, focusing on the white-eye taxa of the island of Lifou near New Caledonia (fig. 1). The subspecies of silvereye whose reference genome we sequenced (Zosterops lateralis melanops) is endemic to Lifou, where it co-occurs sympatrically with two other endemics: Zosterops minutus and Zosterops inornatus (Mees 1957). These Lifou white-eyes exhibit the largest body size differences found in Zosterops, with Z. minutus being the smallest, Z. lateralis melanops intermediate in size, and Z. inornatus being one of the largest species of white-eye (Mees 1957; Black 2010). In order to investigate the pattern, common in the genus, of striking interspecific phenotypic divergence evolving in an oceanic island context within short evolutionary time frames, we used RAD markers of multiple individuals of Z. minutus and Z. inornatus. We also compared the demographic dynamics of the three sympatric Lifou species by means of coalescent analyses at the genome level.

F

Plot of the first two axes from a principal component analysis based on morphological traits measured in sympatric Zosterops species from Lifou Island (a) and geographical location of the Loyalty islands and distribution of the four white-eye taxa studied here (b).

Finally, we investigated genomic evolution of white-eyes at an even finer taxonomic and temporal scale by comparing the genomes of two silvereye subspecies. We contrast the genome of the subspecies of silvereye endemic to Lifou to the resequenced genome of a closely related subspecies—Zosterops lateralis nigrescens—that occurs on two nearby islands, Ouvéa and Maré (fig. 1). The two subspecies exhibit striking differences in the color of head plumage, which is green in Z. lateralis nigrescens and black in Z. lateralis melanops (Mees 1957). By comparing genomes below the species level, we aim to characterize the pattern of differentiation occurring in shallow time that may be responsible for the high diversity of allopatric forms in the silvereye.

In summary, the objectives of this study are to present the first reference genome of a species of white-eye (silvereye) and compare it with 1) genomes of other bird families, in order to investigate the genomic underpinning of rapid speciation in white-eyes; 2) RAD sequences of other species within the genus Zosterops, in order to identify patterns of genomic evolution that may be linked to body size differentiation in the three sympatric endemic taxa of Lifou island, and whole genome resequencing for investigating their past demographic dynamics; and 3) the resequenced genome of a different subspecies of silvereye in order to characterize genomic differentiation that may be associated with intraspecific phenotypic differences.

Materials and Methods

Field Collections and Measurements

We obtained blood samples and morphometric data for white-eyes that inhabit the Loyalty Islands Archipelago, 110 km northeast of New Caledonia (fig. 1b). We collected individuals from the three white-eye taxa of Lifou (167°30′E, 20°45′S), all of which are endemic to the island: The Lifou silvereye (Z. lateralis melanops), the small Lifou white-eye (Z. minutus), and the large Lifou white-eye (Z. inornatus). We additionally collected a different subspecies of silvereye (Z. lateralis nigrescens) found in Ouvéa Island (166°30′E, 20°30′S) and also in the Loyalty Islands. Birds were caught using mist nets during October–December 2011; blood samples were taken via venipuncture from the subbrachial wing vein of each bird and stored either in 95% ethanol, Queen’s lysis buffer, or RNAlater. To avoid resampling, birds were ringed using colored rings. A total of 176 Z. minutus, 87 Z. lateralis melanops, 64 Z. inornatus, and 6 Z. lateralis nigrescens were caught. For each bird, we took the following morphological measurements: Culmen length, depth, and width (to anterior and posterior of nostril); length of the tarso-metatarsus; flattened wing chord length; maximum tail length; head to bill tip length; and body mass. In order to summarize variation in morphology, we performed a multivariate analysis using R.

Genomic DNA Extraction and Genome Sequencing

Genomic DNA was extracted from blood using Qiagen DNeasy Blood and Tissue kits (Qiagen, Inc., USA), from one individual of the silvereye Z. lateralis melanops. The quality and quantity of genomic DNA were assayed using gel electrophoresis and fluoroscence measurements (Qubit). Genome sequencing was performed at BGI, Hong Kong, using an Illumina HiSeq 2000 sequencer. Five paired-end fragment libraries were generated with small insert sizes of 2 × 170, 1 × 500, and 2 × 800 bp plus four longer jumping libraries with insert sizes of 3 × 2,000 and 1 × 5,000 bp. Sequence reads were filtered for quality control as follows. We discarded reads with more than 10% of ambiguous bases or poly-As in the short libraries (170, 500, and 800 bp), and 20% in the long libraries (2 and 5 kb). Reads with more than 40% of low-quality bases (Illumina Q-value < 7) for the short libraries and 30% of low-quality bases for the long libraries were discarded. Reads with adapter contamination (i.e., reads with more than 10 bp aligned to the adapter sequence with no more than 3 bp mismatch allowed) were discarded. Fragment library reads where paired-end reads overlapped by more than or equal to 10 bp (10% mismatch allowed) were discarded. Finally, identical paired-end reads, likely resulting from polymerase chain reaction (PCR) duplicates, were discarded.

Genome Assembly and Quality Check

Genome assembly was performed using ALLPATHS-LG (Broad Institute of MIT and Harvard), a graph-based assembler that has been shown to produce high-quality assemblies for medium- to large-sized genomes (Gnerre et al. 2011). We used version 44.837 with the default parameters, except that the Haploidify option was set to “true,” as recommended for diploid data sets in order to deal with high heterozygosity (switching this option to “false” was found to produce a lower quality assembly). Sequencing depth was obtained by mapping the cleaned fragment reads back to the genome assembly using SOAP2 (Li, Yu, et al. 2009).

A k-mer analysis approach was used to calculate genome heterozygosity and identify repetitive regions. In brief, this method utilizes short sequences of varying length (kbp, 17 bp in our case) from the cleaned short insert size sequence reads. Using a 1-bp sliding window the frequency distribution of each k-mer should be normally distributed, with deviations from this a result of heterozygosity, sequencing error, and repeat regions in the genome.

The genome Guanine-Cytosine (GC) content was also calculated to check for possible sequencing bias of potential contamination.

Genome Annotation

Repetitive sequences in the Zosterops genome were identified and masked with RepeatMasker 4.0.5 (Smit et al. 2010) and RepBase (Jurka et al. 2005) using repeat regions from the chicken genome as a reference. The MAKER 2.10 pipeline (Cantarel et al. 2008) was used for genome annotation with a combined strategy of homology-based and ab initio gene prediction approaches. For the homology-based approach, we aligned protein sequences from chicken, turkey, and zebra finch to the translated nucleotide sequences from Zosterops using TBLASTN (Kent 2002) and Exonerate (Slater and Birney 2005). For the ab initio approach, the SNAP software (Korf 2004) was used to generate a consensus gene set with evidence-based quality values. Only genes with Annotation Edit Distance (AED) score less than 1 were extracted for further analysis.

The final list of genes was then annotated using BLASTP (Camacho et al. 2009) against the NCBI (National Center for Biotechnology Information) nonredundant protein database with matches restricted to the four avian species used for comparative genomics (chicken, turkey, zebra finch, peregrine falcon) with an e-value threshold of 1 × 10−5. BLAST2GO (Conesa et al. 2005) was then used to extract gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes annotations for our genes.

Orthologous Identification

To identify orthologous sequences between avian species, we downloaded the available coding sequences (CDS) of chicken (GCA_000002315.2), turkey (GCA_000146605.1), zebra finch (taeGut3.2.4), and peregrine falcon (http://www.cardiff.ac.uk/people/view/81128-bruford-flsw-mike, last accessed September 2015). We performed a reciprocal best-hit BLAST (BLASTN) using an e-value threshold of 1 × 10−6 implemented in BLAST + version 2.2.28 (Camacho et al. 2009).

Positive Selection and Gene Ontology Enrichment Analysis

To assess the selective pressures acting on the orthologous genes between Zosterops and the other four avian species, we compared the nonsynonymous with synonymous substitution ratio (dn/ds ratio) using likelihood methods implemented in PAML v.4.7 (Yang 2007). Accurate estimates of dn/ds ratio requires high-confidence sequence alignments, as misalignment errors can potentially result in the false positive detection of positive selection (Jordan and Goldman 2012). To generate such alignments, the genes were first translated into protein sequences and then aligned using four different assemblers (MAAFT, MUSCLE, KALIGN, and T_COFFEE) with M-COFFEE (Wallace et al. 2006), part of the T-COFFEE package v.11.0 (Notredame et al. 2000). Consensus alignments from the four different methods were trimmed so that only amino acid residues aligned in identical positions by all assemblers were retained (quality = 9). The trimmed protein alignments were subsequently reverse translated into nucleotide alignments using the original sequences. The heuristic algorithm of MAXALIGN (Gouveia-Oliveira et al. 2007) was then used to identify gap-rich sequences, which may represent paralogous gene sequences. In such cases, these alignments were discarded. To estimate dn/ds ratio for each avian branch in the phylogeny, we ran the “free-ratio” model for a concatenate alignment of the retained genes. Branch-site tests were used to detect episodic positive selection acting on specified “foreground” branches compared with the “background” avian phylogeny (Zhang et al. 2005). Each of the five terminal branches of the phylogeny was individually assigned as foreground branches for the alternative and null models. To ensure convergence of the likelihood methods, all CODEML analyses were ran independently three times with the lowest likelihood value retained (Yang and Dos Reis 2011). Significant positive selection for the site and branch-site models was inferred using log-likelihood ratio tests between the alternative and null models corrected for multiple testing (false discovery rate [FDR] < 0.05; Benjamini and Hochberg 1995). For genes with evidence of positive selection, Bayes Empirical Bayes method was used to assess the significance of the dn/ds ratio at each codon position (posterior probability > 95%; Yang 2007).

The genes potentially evolving under positive selection in the silvereye were used for an enrichment analysis, using the hypergeometric test implemented in GOstat (Falcon and Gentleman 2007).

Genome-Based Phylogeny Reconstruction

We used the 4-fold degenerated sites from genes that did not show evidence of positive selection to build a phylogenetic tree. We used JModelTest 2.1 (Darriba et al. 2012) and identified that GTR + I was the best-fit model based on the Akaike Information Criterion (AIC). Using this substitution model, a maximum-likelihood tree was built using RAxML version 8 (Stamatakis 2014). To date the tree and calculate substitution rates, we used the software MCMCtree from PAML 4.7 (Yang 2007). Because the GTR model is not available in MCMCtree, we used HKY instead. We ran 200,000 iterations (sampling every 2 iterations until the program reached 100,000 samples) after a burn-in of 50,000. We used the age range of the fossil of Vegavis, an extinct bird related to ducks and geese, to set a uniform prior for the root of our phylogenetic tree, that is, between 66 and 86.5 Ma (Benton and Donoghue 2007).

Gene Gains and Losses

The perl script treefam_scan.pl implemented in TreeFam (Li et al. 2006) was used to group protein sequences into gene families. A birth/death process was then modeled in CAFE 3.0 (Han et al. 2013) to estimate gene gains and losses, using an ultrametric tree obtained from MCMCtree. In order to provide an accurate birth/death parameter over the tree (lambda), we removed the gene families that were present in only one species, plus four with large size or variance, as required in CAFE. Possible genome assembly and annotation errors were corrected using the “errormodel” command in the Python script caferror.py.

Detailed Phylogenetic Analyses of Zosterops

To determine whether the three Lifou taxa are the result of a single or multiple colonization events, and to infer when the Loyalty islands were colonized by Zosterops, we constructed a dated phylogenetic tree of all species of Zosterops from the southwest Pacific region surrounding New Caledonia, including taxa that went extinct in historical times. For this we downloaded mitochondrial DNA data (cytochrome B, ND2, ND3, and ATPase) from GenBank (Phillimore et al. 2008; Black 2010; Black, Clegg, Phillimore, Owens, unpublished data). GenBank accession numbers are provided in supplementary table S1, Supplementary Material online. We also sequenced the ND3 gene for some of our own samples of Loyalty island Zosterops, using the PCR protocol described in Black (2010). Phylogeny reconstruction and divergence dating were conducted in BEAST v1.8.0 (Drummond et al. 2012) using the best model of molecular evolution as selected using the AIC in JModelTest 2.1 (GTR + I + G, Darriba et al. 2012). We performed divergence dating using an estimated rate of molecular evolution derived from mitochondrial cytochrome B gene sequences in Passeriformes of 2.07%/Myr (Weir and Schluter 2008). We applied a Bayesian relaxed uncorrelated clock model and ran 4 independent chains for 30 million generations with a birth–death tree prior.

Whole Genome Resequencing

DNA was extracted from one individual for each of the following three taxa: Z. minutus, Z. inornatus, and Z. lateralis nigrescens (Ouvéa Island) and used for library preparation. One paired-end library was constructed with insert size of 500 bp and sequenced on Illumina HiSeq at BGI, Hong Kong. Sequence reads were quality filtered as previously described for the reference genome of Z. lateralis melanops.

RAD Sequencing

We genotyped 2 pools of 10 individuals of Z. minutus and Z. inornatus using RAD sequencing to determine the genomic basis of body size variation. Single-end (100 bp) RADpool libraries containing an equal amount of DNA for each sample per species, prepared with SbfI as restriction enzyme, were processed at Floragenex (University of Oregon High Throughput Sequencing Facility). Reads retained after quality control were aligned to the Z. lateralis melanops genome using Bowtie2 (Langmead and Salzberg 2012) to obtain SAM format alignments. The “view” command in SAMTOOLS 0.1.18 (Li, Handsaker, et al. 2009) was used to convert alignments from SAM to BAM formats and remove ambiguously mapped reads. SAMTOOLS “mpileup” command was then used for generating the input file for POPOOLATION2 (Kofler et al. 2011). POPOOLATION2 was applied for calculating allele frequency differences between Z. inornatus and Z. minutus. For this we set the minimum count of minor allele to 6 and minimum and maximum coverage to 50 and 200, respectively. Fst-sliding.pl perl script was used to identify single nucleotide polymorphisms (SNPs) highly differentiated between species with a sliding window of one. The effect of low coverage on Fst estimates was inspected with the Fisher’s exact test (Fisher 1922) using the fisher-test.pl script from POPOOLATION2. P values were corrected for multiple testing. The significant highly divergent SNPs (Fst > 0.7, FDR < 0.05) were mapped on Z. lateralis melanops genome using the perl script SNPdat_v1.0.5.pl (Doran and Creevey 2013).

Demographic Reconstruction

Historical inferences of effective population size (Ne) for Zosterops taxa were inferred using the pairwise sequentially Markovian coalescent model (PSMC; Li and Durbin 2011). We used SAMTOOLS 0.1.18 (Li, Handsaker, et al. 2009) to generate diploid consensus sequences for each individual with the parameter -C50 to reduce the effect of reads with excessive mismatches. Before running PSMC we converted the diploid sequences into the required format using the perl script fq2psmc.pl. The PSMC model was run for the three Lifou Zosterops taxa with default parameters, 0.5 years as generation time (Moyle et al. 2009) and our previously estimated mutation rate. Reliability of our estimates in Ne fluctuation was evaluated using a bootstrap test run 100 times for each individual using the 100 simulated data sets. For this the genomes were split into small fragments and then randomly sampled, with replacement, using the script splitfa.pl.

Plumage Divergence in Zosterops lateralis

We compared the reference genome of Z. lateralis melanops (black head plumage) with genome resequencing of the closely related subspecies Z. lateralis nigrescens (green head plumage). Reads retained after quality control were synchronized with the reference genome using mpileup command in SAMTOOLS and variants were called with “bcftools.” SNPs obtained were then annotated on Z. lateralis melanops genome using the perl script SNPdat_v1.0.5.pl (Doran and Creevey 2013).

Results and Discussion

The Silvereye Genome

High-quality bird genomic data are rapidly accumulating (Jarvis et al. 2014; Zhang et al. 2014), allowing researchers to address a variety of questions regarding avian evolution that were problematic to tackle prior to the genomic revolution. By presenting the first reference genome for a species of white-eye, we aim to provide a valuable genomic resource for complementing and stimulating new research on various topics in avian biology. We generated 108 Gb of raw Illumina HiSeq high-throughput DNA sequencing data (supplementary table S2, Supplementary Material online), representing approximately 80X coverage of the estimated genome size for a congeneric species (Z. pallidus, 1.35 Gb; Wright et al. 2014). After filtering raw reads for quality control, genome assembly resulted in 80,536 contigs and 2,951 scaffolds with a total length of 1.036 Gb, suggesting a coverage of 64X, and with an N50 of 3.5 Mb and 22.9 Kb for scaffolds and contigs, respectively (supplementary fig. S1 and table S3, Supplementary Material online). The good quality of the genome was confirmed by k-mer and GC content analyses (supplementary figs. S2 and S3, Supplementary Material online). These genomic data have been deposited to GenBank/EBI under the accession LAII00000000.

Comparisons with Other Bird Families

We annotated 20,247 genes in the silvereye genome using both homology assessment with avian species and de novo predictions (Cantarel et al. 2008), resulting in a slightly higher number of genes than other birds (table 1). The proportion of repetitive elements was similar to that observed in other avian species and accounted for 5.61% of the genome (supplementary table S4, Supplementary Material online). A reciprocal best-hit BLAST (e-value < 10−6) approach identified 6,730 potentially orthologous genes among CDS of the five species analyzed. Further analyses using multiple assemblers and removing putative paralogs resulted in high-quality consensus alignments of 3,682 orthologous genes. From the 1,511 genes that did not show evidence of positive selection (see below), we obtained 254,348 4-fold degenerate sites that were used to build a phylogenetic tree and calculate substitution rates. The tree, calibrated with a fossil for the zebra finch–chicken split (Benton and Donoghue 2007), showed divergence times between species in concordance with previous phylogenetic analyses (Zhan et al. 2013). Furthermore, our tree shows that zebra finch and silvereye diverged approximately 23 Ma (34-12 My, 95% confidence interval [CI]), more recently than previously identified (44.7 My, fig. 2 and supplementary fig. S4, Supplementary Material online; Jetz et al. 2012). The degenerate sites were also used to calculate substitution rates: As expected from our hypotheses, we found that the branch leading to the silvereye had the highest substitution rate (3.16 × 10−9 substitutions/site/year). However, note that because of the relatively deep split between the silvereye and the zebra finch, some of the substitutions that have accumulated along the branch leading to Zosterops may have arisen before this genus split from its sister genus. We also confirmed that the substitution rates were lowest in the falcon (1.76 × 10−9, as reported by Zhan et al. 2013), while turkey, chicken, and zebra finch presented intermediate rates (1.99, 1.88, and 2.25 × 10−9 substitutions/site/year, respectively; supplementary table S5, Supplementary Material online).

Table 1

Summary Statistics of Bird Genomes

	Genome Assembly Size (Gb)	Percentage GC Content (%)	Number of Genes	Percentage of Repetitive Elements (%)	Scaffold Length (N50, Mb)	Platform
Zebra finch	1.196	41.3	18,204	7.10	10.00	Sanger
Chicken	1.031	41.2	16,354	8.56	7.07	Sanger
Peregrine falcon	1.174	41.7	16,263	4.30	3.89	Hiseq 2000
Turkey	1.061	40.5	16,496	7.49	1.50	GAII, 454
Silvereye	1.036	41.1	20,247	5.61	3.50	Hiseq 2000

F

Phylogeny and genome evolution in white-eyes relative to other birds. On the left, genome-based phylogenetic tree showing time scale in million years before present, rate of gene gains (+) and losses (−) per million year (above branch), substitution rate (substitution/site/year × 10−9; below branch, left), and branch-specific dn/ds ratio (below branch, right). On the right, mitochondrial DNA-based phylogenetic tree of the white-eyes of Lifou, Ouvéa, and Maré, showing time scale in million years before present, geographic distribution, body mass (above branch), and phenotypes (photos, O. Hébert).

Gene duplication is considered one of the main genetic sources for evolutionary innovations, having been shown to fuel biological diversification and speciation (Lynch 2000; Brawand et al. 2014). Over 90% of genes for each species were classified into families, with an average of 2.15 loci per gene family (supplementary table S6, Supplementary Material online). Again, as we expected, the branch leading to the silvereye showed the highest frequency of gene duplication, with an average of 0.30 genes gained per family (table 2). This is indicative of a conspicuous high rate of gene duplication in this part of the tree, with 133 families showing significant changes in gene numbers (FDR < 0.05; 127 gains and 6 losses; fig. 2 and table 2). Gene gains were observed in gene families involved in biological processes such as cell adhesion (teneurin, cadherin), cell/cell interaction (ephrin, fibronectin), locomotion (myosin, nesprin), and signal transduction (tyrosine kinase, glutamate receptor; supplementary table S7, Supplementary Material online).

Table 2

Comparison of the Evolution Gene Families in Avian Genomes

	Zebra finch	Silvereye	Falcon	Chicken	Turkey
Number of gene family	7,854	7,885	7,745	7,916	7,603
Average expansion (mean number of genes gained [+] or lost [−] per family)	+0.023	+0.297	−0.188	−0.052	+0.049
Families showing an expansion	1,272	1,785	445	559	1,127
Families remained constant	5,623	5,712	5,894	6,598	616
Families showing a contraction	1,313	712	1,870	1,052	922
Gene gained per million year	55.30	77.61	8.24	22.36	45.08
Gene lost per million year	−57.09	−30.96	−34.63	−42.08	−36.88

To look for possible positive selection acting upon orthologous genes identified across avian genomes, we calculated the ratio of nonsynonymous to synonymous substitutions (dn/ds) on a set of highly reliable genes (3,682, 55% of orthologous genes). Branch-specific dn/ds ratio in the silvereye lineage was above average (0.173), albeit not the highest (fig. 2), but nevertheless indicating accelerated functional evolution in Zosterops. Because the branch leading to Zosterops is relatively long (∼20 Ma), the signature of selection acting at the start of this branch may have been lost by subsequent neutral evolution. The calculations of dn/ds ratios for each individual locus across the phylogeny showed that the branch leading to the silvereye had the highest number of genes potentially evolving under positive selection (990 genes with dn/ds > 1, P < 0.05; 881 genes with FDR < 0.05), of which 525 are unique to this lineage (table 3).

Table 3

Number of Genes Potentially Evolving under Positive Selection (See Methods for Details)

Species	P Value	FDR	Bonferroni	Exclusive
Silvereye	990	881	576	525
Chicken	140	59	40	24
Falcon	437	342	235	183
Turkey	542	489	369	271
Zebra finch	588	466	302	240

Among the genes under positive selection along the branch leading to the silvereye, we conducted GO enrichment analyses, which, after correcting for multiple testing (FDR < 0.20), detected that 11 “biological processes,” 15 “cellular components,” and 27 “molecular functions” were significantly enriched (supplementary table S8, Supplementary Material online). We found that 9 of the 11 enriched biological processes were related to actin filament-based or regulation of actin polymerization/depolymerization processes (FDR < 0.20, supplementary table S8, Supplementary Material online, and list of genes in table 4). Actin, one of the major cytoskeletal components in eukaryotic cells, provides a link to both dispersal abilities and reproductive isolation (Rayment et al. 1993; Pollard and Cooper 2009). First, actin, together with other molecules (i.e., troponin; table 4), plays a key function in muscle contraction, and thereby in flight. High frequency of muscle contractions is necessary to produce enough power to overcome air resistance and enable flight, especially for small-sized birds (Biewener 2011). Second, actin polymerization occurs during sperm capacitation followed by a rapid depolymerization during the acrosome reaction in mammals (Breitbart et al. 2005). Actin polymerization is regulated by actin-related protein 2 (Lee et al. 2013), while destrin or gelsolin is involved in its depolymerization (table 4; Howes 2001; Finkelstein et al. 2013). Positive selection on genes related to actin poly/depolymerization supports the idea that reproductive proteins evolve rapidly during the establishment of reproductive barriers that lead to speciation (Swanson and Vacquier 2002). Hence, we suggest a plausible mechanism by which actin-related genes may facilitate island colonization (through muscle enhancement) while erecting reproductive barriers (through fertilization incompatibility), thus providing a potential solution to the paradox of the great speciator. Although speculative, these results indicate that pleiotropy and evolutionary flexibility might hold the key to rapid speciation in white-eyes.

Table 4

List of Genes in the Silvereye that are Both Potentially Evolving Under Positive Selection and Significantly Enriched for Biological Processes

Transcript ID	Gene Description
ENSGALT00000039770	Spectrin, alpha, nonerythrocytic 1 (alpha-fodrin, SPTAN1)*
ENSGALT00000002199	Gelsolin (GSN)*
ENSGALT00000014276	Actin-related protein 2 (ACTR2)*
ENSGALT00000018647	Villin 1 (VIL1)*
ENSGALT00000014113	Destrin (actin depolymerizing factor, DSTN)*
ENSGALT00000015293	Capping protein (actin filament) muscle Z-line, alpha 2 (CAPZA2)*
ENSGALT00000045497	ADP-ribosylation factor 6 (ARF6)*
ENSGALT00000004183	Adapter molecule crk*
ENSGALT00000044680	Neurotrophin 3 (NTF3)*
ENSGALT00000002214	Troponin C type 1 (TNNC1)*
ENSGALT00000009894	Formiminotransferase cyclodeaminase (FTCD)
ENSGALT00000039415	Lactate dehydrogenase A (LDHA)
ENSGALT00000008125	Prostaglandin-endoperoxide synthase 2 (PTGS2)
ENSGALT00000015194	Degenerative spermatocyte homolog 1 (DEGS1)
ENSGALT00000019819	ELOVL fatty acid elongase 6 (ELOVL6)
ENSGALT00000002022	Basigin (BSG)

Note.—Asterisk indicates genes involved in actin filament-based or regulation of actin polymerization/depolymerization processes.

Our comparison of Zosterops with other bird families consistently supports the view that white-eye genomes likely present unusual evolutionary flexibility relative to the other avian lineages included in our analyses, a feature that may have provided the necessary source material for rapid speciation and phenotypic differentiation. Although we acknowledge that including a wider range of avian genomes will be required in order to confirm this hypothesis, it is clear that our results are consistent with “high evolvability” in Zosterops. Future work should not only aim at a broader scope in the sampling of avian lineages, but also at clearly pinpointing the key genes involved, with a specific focus on actin-related genes that we found to be of particular relevance.

Comparisons within Zosterops

For our analyses of genome differentiation between different species of Zosterops, we focused on the white-eyes of Lifou island, where three species occur in sympatry, whereas most islands tend to support only one or two species (Gill 1971). Our morphometric analyses revealed that the three Lifou taxa exhibit significant, striking differences in body size: Z. minutus has an average body mass of 9 g (7–11 g, 95% confidence interval), the silvereye 14 g (13–16 g), and Z. inornatus 22 g (19–25 g) (supplementary table S9, Supplementary Material online). These are exceptional in terms of body size differences, including the two extremes in the genus: Z. minutus, being the smallest, and Z. inornatus being the largest, both are found in sympatry (Mees 1957; Black 2010). The phylogenetic reconstruction including complete sampling of the entire diversity of white-eyes in the southwest Pacific region (Black 2010) showed that Z. lateralis melanops, Z. minutus, and Z. inornatus have independently colonized Lifou, therefore excluding sympatric speciation on Lifou (supplementary fig. S5, Supplementary Material online). Calibrating the tree with the estimated rate of molecular evolution for mitochondrial sequences of passeriformes, we found that Z. inornatus was the first colonizer of the island (1.50 Ma, 0.95–2.07 Ma 95% highest posterior density interval), followed by Z. minutus 1.15 Ma (0.74–1.58) and Z. lateralis (less than 0.5 Ma; supplementary fig. S5, Supplementary Material online). This is consistent with the idea that avian body size tends to change following island colonization (Clegg and Owens 2002; Frentiu et al. 2007).

In order to gain insights into a typical pattern of genomic differentiation in the white-eye radiation, we searched for the genomic underpinning of body size differences in the Lifou taxa. Our phylogenetic analysis revealed that Z. minutus and Z. inornatus diverged about 2.06 Ma (1.67–2.47), suggesting that the significant divergence in body size between these species has evolved recently. This clade (i.e., including Z. minutus and Z. inornatus) split from the Z. lateralis clade 2.77 Ma (2.41–3.14) (supplementary fig. S5, Supplementary Material online). We genotyped pools of ten individuals of Z. minutus and Z. inornatus using RAD sequencing with the aim of identifying genomic regions with high level of divergence between the two species. We detected 5,245 SNPs with high divergence (Fst > 0.7) and significant differences in allele frequency (Fisher’s exact test, FDR < 0.05). We did not do a GO enrichment analysis, but instead mapped these SNPs onto the silvereye genome; we found that 185 of them are within the exons of 162 genes. Despite RAD sequences representing a small portion of the genome, we nevertheless identified 5 genes among these 162 that are potentially good candidates underlying body size divergence, being involved in signal transduction (“inositol 1,4,5-trisphosphate receptor type 2 and 3”, both with a key role in metabolism and growth) and structural (“collagen 6A3” and “17A1”), as well as tissue (“keratin type I”) development (supplementary table S10, Supplementary Material online).

We also conducted a genome-wide coalescence analysis in order to estimate fluctuation in effective population sizes of the three Lifou taxa, which may be indicative of past interspecific interactions in a context of sympatry. Our reconstruction revealed that these taxa have experienced contrasting demographic histories. In particular, during the last 100,000 years, Z. inornatus and Z. minutus underwent population expansions while the silvereye saw a reduction (fig. 3). About 25,000 years ago the effective population size of Z. inornatus was reduced by half, while that of Z. minutus increased and that of the silvereye shrank. Although the low-depth coverage of the resequenced genomes could have an effect on the estimated effective population sizes, the demographic dynamic patterns should not be affected by coverage (Miller et al. 2012). The contrasting demographic trends observed for the three sympatric species of Lifou could be due to high interspecific competition in Zosterops. The reconstructed demographic trajectories might be compatible with a dynamic scenario of successive island colonizations followed by competitive exclusion, which potentially set the stage for the striking divergence in body size between the three sympatric species that we observe today.

F

Demographic history of the white-eyes of Lifou island based on coalescence analyses of their genomes. The bold lines show population sizes through time, with 100 bootstrap analyses indicated by thin lines. Orange: Zosterops lateralis melanops; red: Zosterops inornatus; green: Zosterops minutus.

Our analyses of white-eye genomes of different species from Lifou island reveal that high levels of interspecific differentiation have accumulated in short evolutionary time frames and that multiple genomic regions appear to have been linked to body size divergence. In addition, the detected signal of strong contrasting past demographic fluctuations supports a scenario of negative competitive interactions within a single island, thus providing insights into why only fewer than two species of Zosterops tend to co-occur in a single island (Gill 1971). Denser taxonomic sampling of the white-eye clade that includes Z. minutus and Z. inornatus, combined with higher depth of coverage of genomic sequence data for these taxa may reveal precisely what regions are linked to the remarkable body size differentiation between these closely related bird species.

Within-Species Comparisons

The two subspecies of silvereye (Z. lateralis melanops and Z. lateralis nigrescens) diverged less than half a million year ago (supplementary fig. S5, Supplementary Material online) and exhibit no differences in body size (fig. 1a), but striking differences in head color (supplementary fig. S6, Supplementary Material online). The alignment of the resequenced Z. lateralis nigrescens genome to the Z. lateralis melanops reference genome revealed 3,700,336 SNPs, 124,105 of which mapped to exonic gene regions and 25,374 of them resulted in nonsynonymous substitutions. The 6,737 Zosterops genes containing at least one nonsynonymous substitution were blasted against the chicken CDS gene set. Among the BLAST best-hits (e-value < 10−6), we found that 12 genes were related to bird feather color through the melanin or carotenoid synthesis pathway (table 5); none of them resulted to be under positive selection. Of particular interest was MIFT, known to affect melanin synthesis and the color of bird plumage (Wang et al. 2014), and which was truncated by a stop codon in Z. lateralis nigrescens. Although more detailed analyses would be necessary to obtain exhaustive conclusions, our analysis highlights some genomic regions that can provide a good primer for further tests of the genetic basis of green and black-headed plumage color in the silvereye, which may in the future help explain more widely patterns of melanism in white-eyes and other birds.

Table 5

List of Plumage-Linked Genes in Various Bird Orders with Indication of Numbers of Nonsynonymous Substitutions in Zosterops lateralis nigrescens When Compared with Zosterops lateralis melanops

Gene	Number of Nonsynonymous Substitutions	Reference	Pathway
Tyrosinase	9	Xu et al. (2013), Wang et al. (2014)	Melanin
Tyrosinase-related protein 1	1	Xu et al. (2013)	Melanin
Microphthalmia-associated transcription factor	1	Wang et al. (2014)	Melanin
Cluster of differentiation 36	3	Walsh et al. (2012)	Carotenoids
Beta-carotene 15,15′-monooxygenase 1	1	Walsh et al. (2012)	Carotenoids
StAR-related lipid transfer	39	Walsh et al. (2012)	Carotenoids
Endothelin receptor B	3	Ekblom et al. (2012)	Melanin
Endothelin receptor B2	7	Kinoshita et al. (2014)	Melanin
Agouti-related peptide	1	Li et al. (2011), Yoshihara et al. (2012)	Melanin
G-protein coupled receptor 143	2	Shi et al. (2012)	Melanin
Paired box protein-3	1	Li et al. (2011)	Melanin
Cyclin-dependent kinase inhibitor	5	Hellström et al. (2010)	Melanin

In conclusion, we have conducted a set of genomic analyses at different taxonomic scales in order to gain insights into the genetic basis of fast diversification in white-eyes. Given that Zosterops is one of the most rapidly diversifying lineages of birds (Moyle et al. 2009; Jetz et al. 2012), we hope that the silvereye genome will help shed light on the drivers of species radiations, as well as on the ability of birds to colonize and adapt to newly invaded environments. Although our study has provided some clues as to how the paradox of the great speciator may be solved, for example, through the rapid evolution of actin-related genes, further analyses including multiple individuals from the many islands colonized by Zosterops would be welcome in the future.

Supplementary Material

Supplementary figures S1–S6 and Tables S1–S10 are available at Genome Biology and Evolution online (http://www.gbe.oxfordjournals.org/).