Journal of Heredity, 2020, 21–32
doi:10.1093/jhered/esz057
Symposium Article
Advance Access publication November 14, 2019
Symposium Article
Adaptive Radiation Genomics of Two 
Ecologically Divergent Hawai‘ian Honeycreepers: 
The ‘akiapōlā‘au and the Hawai‘i ‘amakihi
Michael G. Campana, André Corvelo, Jennifer Shelton,  
Taylor E. Callicrate, Karen L. Bunting, Bridget Riley-Gillis, Frank Wos,  
Justin DeGrazia, Erich D. Jarvis, and Robert C. Fleischer
From the Center for Conservation Genomics, Smithsonian National Zoo and Conservation Biology Institute, 
Washington, DC 20008 (Campana, Callicrate, and Fleischer); New York Genome Center, New York, NY 10013 (Corvelo, 
Shelton, Bunting, Riley-Gillis, Wos, and DeGrazia); Species Conservation Toolkit Initiative, Chicago Zoological 
Society, Brookfield, IL 60513 (Callicrate); The Rockefeller University, New York, NY 10065 and Howard Hughes 
Medical Institute, Chevy Chase, MD 20815 (Jarvis).
Address correspondence to Michael G. Campana at the address above, or e-mail: campanam@si.edu.
Received March 4, 2019; First decision May 9, 2019; Accepted September 30, 2019.
Corresponding Editor: Rosemary Gillespie
Abstract
The Hawai‘ian honeycreepers (drepanids) are a classic example of adaptive radiation: they adapted 
to a variety of novel dietary niches, evolving a wide range of bill morphologies. Here we investigated 
genomic diversity, demographic history, and genes involved in bill morphology phenotypes in 2 
honeycreepers: the ‘akiapōlā‘au (Hemignathus wilsoni) and the Hawai‘i ‘amakihi (Chlorodrepanis 
virens). The ‘akiapōlā‘au is an endangered island endemic, filling the “woodpecker” niche by 
using a unique bill morphology, while the Hawai‘i ‘amakihi is a dietary generalist common on 
the islands of Hawai‘i and Maui. We de novo sequenced the ‘akiapōlā‘au genome and compared 
it to the previously sequenced ‘amakihi genome. The ‘akiapōlā‘au is far less heterozygous and 
has a smaller effective population size than the ‘amakihi, which matches expectations due to its 
smaller census population and restricted ecological niche. Our investigation revealed genomic 
islands of divergence, which may be involved in the honeycreeper radiation. Within these islands 
of divergence, we identified candidate genes (including DLK1, FOXB1, KIF6, MAML3, PHF20, RBP1, 
and TIMM17A) that may play a role in honeycreeper adaptations. The gene DLK1, previously shown 
to influence Darwin’s finch bill size, may be related to honeycreeper bill morphology evolution, 
while the functions of the other candidates remain unknown.
Subject areas: Genomics and gene mapping, Conservation genetics and biodiversity
Keywords: bill morphology, Chlorodrepanis virens, demography, Hemignathus wilsoni, islands of divergence, natural selection
Adaptive radiation is a process that involves branching speciation The genomics of this process has only recently begun to be studied 
with differential adaptation (Futuyma 1986). It usually involves mor- (Berner and Salzburger 2015; Wolf and Ellegren 2017; Malinsky 
phological change that allows species to occupy alternative niches, et al. 2018; Salzburger 2018), with some key examples deriving from 
and derived species may become more ecologically specialized. studies of groups such as African lake cichlids (Brawand et al. 2014; 
Published by Oxford University Press on behalf of The American Genetic Association 2019. This work is written by (a) US Government 21
employee(s) and is in the public domain in the US.
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
22 Journal of Heredity, 2020, Vol. 111, No. 1
Henning and Meyer 2014; Malinsky et al. 2018; Salzburger 2018), and H. affinis: Fleischer 2009; Fleischer RC, Campana MG, et al., 
Heliconius butterflies (Joron et al. 2011; Reed et al. 2011; Supple unpublished data), the ‘akiapōlā‘au has a unique bill morphology—
et al. 2013; Kozak et al. 2018), and Anolis lizards (Tollis et al. 2018). the nukupu‘us had longer bills and a curved lower mandible, while 
Of particular interest is discovering the genomic architectures that the ‘akiapōlā‘au has a straight and stout lower mandible (Pratt 
enable adaptive radiation to proceed, and the number and types of 2005). The nukupu‘us apparently became extinct by the late 1990s 
genes and genetic changes involved in the processes of speciation (although accuracy of sightings post-1900 is debated: Elphick et al. 
and adaptation (e.g., Salzburger 2018; Tollis et al. 2018). In birds, 2010; Roberts et  al. 2009). The ‘akiapōlā‘au is itself endangered 
the latter aspect has gained insight via the study of the radiation of (International Union for the Conservation of Nature [IUCN] Red 
Darwin’s finches (Geospizinae), through both gene expression ana- List categories: B1ab(i,ii,iii,iv,v); C2a(ii)) (IUCN 2019), necessitating 
lyses (Abzhanov 2004, 2006), genome scans and genome-wide as- conservation action to prevent its extinction (e.g., Gorresen et  al. 
sociation studies (GWAS) (Lamichhaney et  al. 2015, 2016). These 2007).
studies identified genes (ALX1 and BMP4) and the calmodulin Unfortunately, the Hawai‘ian honeycreepers in general are an 
pathway as being critical to the adaptation of Darwin’s finches’ bill extremely threatened clade: of the 17 confirmed extant species (18 
morphologies to varying environments and diets. Unfortunately, few if Maui ‘amakihi [Chlorodrepanis virens wilsoni] is considered a 
other avian radiations have been analyzed using such methods to separate species rather than subspecies), the IUCN rates 6 as critic-
date. ally endangered, 4 as endangered, and 5 as vulnerable. Another 22 
The Hawai‘ian honeycreepers (drepanids) are a spectacular ex- species (and the Lānaʻi ʻalauahio subspecies [Paroreomyza montana 
ample of rapid adaptive radiation in a species-depauperate island montana]) either went extinct or likely went extinct (no recent con-
system. Currently classified in the Carduelinae subfamily of the firmed sightings) since the arrival of Europeans to the Hawai‘ian 
Fringillidae (finch) family (Fleischer et al. 2001; James 2004; Lerner islands. At least 17 additional extinct species are known only from 
et  al. 2011; Zuccon et  al. 2012), the more than 50 honeycreeper subfossil material (James and Olson 1991). In addition to the grave 
species arose from a single colonization of the remote Hawai‘ian is- conservation concern for these emblematic Hawai‘ian forest birds, 
lands by rosefinch-like (Carpodacus spp.) ancestors between 5.2 and the endangerment and extinction of most members of the clade hin-
8.4 million years ago (Mya; Lerner et al. 2011). In addition to the ders our ability to research most of these taxa due to the limited 
typical “finch” granivorous niche (e.g., palila [Loxioides bailleui]), availability of samples.
honeycreepers fill nectarivorous (e.g., ‘apapane [Himatione Callicrate et al. (2014) published a reference-assisted de novo as-
sanguinea]), generalist (e.g., Hawai‘i ‘amakihi [Chlorodrepanis sembly of the Hawai‘i ‘amakihi genome. The Hawai‘i ‘amakihi is 
virens]), arthropodivorous (e.g., ‘akikiki [Oreomystis bairdi]), and a dietary generalist honeycreeper species endemic to the islands of 
frugivorous (e.g., ‘o‘u [Psittirostra psittacea]) roles. To fill these eco- Hawai‘i and Maui. Each island has a unique subspecies, and the 2 
logical niches, honeycreepers evolved numerous bill morphologies, subspecies are genetically differentiated (Tarr and Fleischer 1993; 
including some relatively unique forms of unknown function in Fleischer et al. 1998; Lerner et al. 2011). Rated as Least Concern 
some of the extinct species (James 2004), the genetics behind which by the IUCN, the Hawai‘i ‘amakihi is common throughout its range 
are currently unknown. (Pratt 2005), even occupying lower elevations on Hawaii Island des-
One of the most extreme examples of bill morphological adap- pite the presence of invasive avian malaria (Plasmodium relictum 
tation is the ‘akiapōlā‘au (Hemignathus wilsoni; Figure 1), the last strain GRW4; Beadell et al. 2006)—lethal to most honeycreepers—
confirmed extant species in its genus. The ‘akiapōlā‘au is endemic to in these regions (Foster et  al. 2007; Atkinson et  al. 2013). It has 
the dry and montane moist forests of the youngest island of Hawai‘i evolved a relatively small and thin bill and feeds largely on a mixed 
(Ralph and Fancy 1996). An ecological specialist, the ‘akiapōlā‘au diet of invertebrates and nectar (Pratt 2005; Figure 1).
fills the “woodpecker” niche, using its shorter, straight mandible to Currently, we can only compare the ‘amakihi genome with dis-
peck bark and its longer, curved maxilla to extract insects for con- tantly related nondrepanid passerine genomes, severely limiting 
sumption. The species prefers to forage on koa (Acacia koa) trees, our ability to detect traits related to the honeycreeper adaptive 
but will also utilize ʻōhiʻa lehua (Metrosidros polymorpha) and naio radiation. Additional genomic resources (including genomes, tran-
(Myoporum sandwicense) (Ralph and Fancy 1996). It has small scriptomes, exomes, genome-wide marker panels, etc.) are critical 
census population sizes (currently ~1000 individuals: Gorresen et al. to understand the honeycreeper radiation and evolution of bill 
2007) due to its limited niche and slow population growth rates morphologies. These resources would permit the interrogation of the 
(Ralph and Fancy 1996). Even compared to its closest known re- complete genome (including regulatory regions, putative noncoding 
latives, the extinct nukupu‘us (Hemignathus hanapepe, H. lucidus, sequences, and genes of unknown function) for signs of selection 
Figure 1. Comparison of bill morphologies between the ‘akiapōlā‘au (Hemignathus wilsoni) and the Hawai‘i ‘amakihi (Chlorodrepanis virens). The specialist 
‘akiapōlā‘au uses its straight mandible to peck bark and its curved maxilla to extract insects. The ‘amakihi has a relatively small and thin bill and feeds largely on 
a mixed diet of invertebrates and nectar. Photographs by Jack Jeffrey (used with permission).
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
Journal of Heredity, 2020, Vol. 111, No. 1 23
(e.g., greatly increased or decreased genomic heterozygosity and di- The genome was assembled using ABySS 1.9.0 (Simpson et al. 
vergence), thereby greatly increasing detection power for causative 2009) with a K-mer length of 71. Gap closure was performed using 
traits over candidate gene approaches. For instance, only with the ABySS 1.9.0 Sealer under default settings except that we used a larger 
power of transcriptomics, genome scans, and GWAS were genes re- bloom filter size (2400 MB) and maximum number of branches in 
sponsible for the variation in bill morphologies in Darwin’s finches the de Bruijn graph traversal (3000). We assembled the genome 
identified (Abzhanov 2004, 2006; Lamichhaney et al. 2015, 2016). using a range of K-mer sizes (61, 71, 75, 79, 83, 87) and chose the 
Similarly, the 4 candidate regions that likely explain plumage dif- assembly that produced the largest scaffold N50 values.
ferences between golden-winged (Vermivora chrysoptera) and blue- Genome assembly statistics were calculated using the 
winged warblers (V. cyanoptera) required analysis of these species’ assemblathon_stats.pl script from the Assemblathon 2 competition 
complete genomes (Toews et al. 2016). (Bradnam et  al. 2013). Genic completeness was calculated using 
Therefore, to explore genomic effects of niche specialization, iden- BUSCO 3.0.2 (option --long) with the aves_odb9 gene set (totaling 
tify candidate genes and regions to explain the rapid radiation of 4915 orthologs) and the AUGUSTUS “chicken” model as a starting 
honeycreeper bill morphology, and inform honeycreeper conservation data set (Stanke et  al. 2006a; Waterhouse et  al. 2017). For com-
strategies, we de novo sequenced the genome of the ‘akiapōlā‘au. The parison, we also ran BUSCO on the ‘amakihi genome (both the 
recent divergence of the Hawai‘i ‘amakihi from the ‘akiapōlā‘au (~3.0 Callicrate assembly and the ‘akiapōlā‘au-realigned sequence; see 
Mya: Lerner et al. 2011) facilitates genomic alignment. In addition, below).
their co-occurrence on the island of Hawaii, which provides a max-
imum age of about 1 million years (Lipman and Calvert 2013), and Variant Calling and Consensus Sequence 
their differing degrees of morphological and ecological specialization, Generation
makes the comparison of the genomes of the 2 species potentially im-
The Callicrate ‘amakihi genome sequence is not directly compar-
portant for determining patterns and rates of genome evolution and 
able to our de novo ‘akiapōlā‘au genome due to differing sequencing 
adaptation. We tested the following hypotheses:
technologies and assembly approaches. The ‘amakihi contigs were 
1. Due to its smaller census size and limited niche, the ‘akiapōlā‘au generated from a combination of ~2× Roche/454 and ~19× Illumina 
will have lower heterozygosity, a smaller effective population GAII sequences (Callicrate et al. 2014). The authors generated ~60× 
size (N ) and a differing demographic history from the more nu- Illumina coverage, but downsampled their data for contig assembly e
merous, dietary generalist ‘amakihi. due to both computer memory limits and the need to maintain the 
2. The fully de novo ‘akiapōlā‘au genome sequence will facili- Roche/454 long reads in the assembly. The contigs were then oriented 
tate more accurate variant calling and genomic reconstructions to the zebra finch (Taeniopygia guttata) genome (Warren et al. 2010) 
within honeycreepers than analyses using divergent reference to generate a chromosome-level assembly. This genomic orientation 
sequences. could impact genomic analyses and interpretations by incorrectly 
3. The ‘akiapōlā‘au and ‘amakihi genomes will be more divergent assuming synteny where there have been genomic rearrangements.
at or near loci responsible for adaptive traits such as bill morph- Therefore, in order to compare the genomes of the ‘akiapōlā‘au 
ology than in the rest of the genome. These divergent regions may and the ‘amakihi genomes, we realigned the short-insert Illumina 
be in previously identified adaptive genes that may be convergent reads from both species to the ‘akiapōlā‘au genome assembly using 
in other species, newly discovered genes of unknown function, BWA-MEM 0.7.12 (Li 2013). We omitted the ‘akiapōlā‘au mate-
putatively noncoding sequences, or regulatory regions. pair and ‘amakihi Roche/454 sequences to minimize batch effects 
due to differing sequence technologies. PCR duplicates were re-
moved using Picard 1.141 MarkDuplicates (Broad Institute 2015). 
Methods Sequence variants were called and consensus sequences were gener-
ated using SAMtools 1.2 mpileup (option -C50) and BCFtools 1.2 
Genome Sequencing and Assembly call (option -c) (Li et al. 2009). We used the older consensus variant 
We sequenced a wild male ‘akiapōlā‘au (SOL748, USGS band caller since downstream PSMC (pairwise sequential Markovian co-
number 1581-74856) captured at Solomon’s Camp, Mauna Loa, alescent) analysis (Li and Durbin 2011) is not compatible with the 
Hawaii, United States on 22 February 2002. Whole blood was newer multiallelic algorithm. Variant calls were converted to FASTQ 
obtained by brachial venipuncture. Genomic DNA was extracted format using vcfutils.pl vcf2fq. We retained single-nucleotide vari-
using a standard phenol-chloroform isolation protocol (Tarr and ants (SNVs) with quality scores of at least 20 for subsequent analysis.
Fleischer 1993). Three Illumina libraries were constructed from the To control for ‘akiapōlā‘au reference bias (defined here as false 
genomic DNA at the New York Genome Center: a standard paired- positive and negative variants derived from inaccurate mapping 
end short-insert (550 base pair [bp] average length) library and 2 to a divergent reference sequence) impacting our downstream esti-
mate-pair libraries with insert lengths of 3–5 kilobase pair (kbp) mates of substitution rates and demography reconstructions, we also 
and 8–10 kbp, respectively. Libraries were 2 × 125 bp paired-end aligned the ‘akiapōlā‘au and ‘amakihi short-insert Illumina reads to 
sequenced on 3 lanes of an Illumina HiSeq 2500. the zebra finch genome (assembly taeGut3.2.4; Warren et al. 2010) 
Polymerase chain reaction (PCR) duplicates were removed using and the Callicrate ‘amakihi genome assembly using BWA-MEM 
an in-house tool at the New York Genome Center. Adapters and low- 0.7.17 (Li et al. 2013). Mate-pair tags were fixed and PCR duplicates 
quality bases were trimmed using Cutadapt 1.8.1 (Martin 2011). For were marked using SAMtools 1.9 (Li et al. 2009) fixmate (option -m) 
the mate-pair libraries, junction adapter identification and clipping and markdup, respectively. Indels were left aligned using the Genome 
were also performed using Cutadapt. Contaminants were removed Analysis Toolkit 4.1.0.0 LeftAlignIndels command (McKenna et al. 
via GEM pre-rel3 by alignment against a spiked-in PhiX sequencing 2010). Variants were called using BCFtools 1.9 mpileup (option -a 
control (Marco-Sola et  al. 2012). Joint error correction was per- DP,AD) and BCFtools call (option -c) (Li et al. 2009). We retained 
formed using Lighter 1.0.6 (Song et al. 2014). SNVs with quality scores of at least 20.
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
24 Journal of Heredity, 2020, Vol. 111, No. 1
Genomic Diversity and Divergence 2 significant figures and used an autosomal genome length of 950 
Because the ‘akiapōlā‘au was male and the ‘amakihi was female, we Mbp for all analyses.
excluded the sex chromosomes from downstream heterozygosity, di- As a first estimate of autosomal substitution rates, we assumed 
vergence, and demographic analyses. Moreover, exclusion of the Z complete lineage sorting within the honeycreeper tree and that all 
improved these estimates since the sex chromosomes are difficult to observed within-species polymorphisms arose after colonization of 
reconstruct accurately due to the presence of repeat elements and Hawai‘i island. Although this assumption could be incorrect (incom-
homologous sequences between the Z and W (Smeds et al. 2015). We plete lineage sorting is common in birds: Jarvis et al. 2014), the pat-
also excluded the nonrecombining mitochondrial genome due to its terns of differentiation (Tarr and Fleischer 1993; Eggert et al. 2009) 
maternal inheritance. For the zebra-finch- and Callicrate-‘amakihi- and extremely low degree of hybridization among honeycreeper 
aligned analyses, we excluded the previously annotated mitogenome taxa (Knowlton et al. 2014) suggest that gene flow from Maui or 
and Z chromosome assemblies (Z and Z_random). For the de novo other older islands was not likely for either species and probably 
‘akiapōlā‘au-aligned analyses, we identified putative mitochondrial ceased once they became differentiated from their sister taxa. Since 
and Z chromosome scaffolds by searching the ‘akiapōlā‘au genome each mutation could occur on either parental chromosome copy 
for matches to the previously sequenced ‘akiapōlā‘au mitogenome (Malinsky et al. 2018), we estimated the within-species substitution 
(GenBank accession KM078774.2; Lerner et al. 2011) and the zebra rates using the following formula:
finch Z chromosome assembly (taeGut 3.2.4; Warren et  al. 2010) 
using megaBLAST (BLASTN 2.7.1+: Camacho et al. 2008). Scaffolds µC = NH/L/TC/2
with matches of bit score 4000 or greater, sequence identity 90% or 
greater, and alignment length 1000 or greater were considered to Here µ C is the genomic single-nucleotide substitution rate, NH is the 
derive from the mitogenome or Z chromosome. We excluded 3 mito- number of heterozygous single-nucleotide polymorphisms within the 
chondrial (1127, 1842, 2102) and 59 Z chromosome (28, 47, 55, 91, branch, L is the sequence length, and TC is the time of colonization. 
101, 105, 107, 111, 120, 122, 131, 133, 141, 144, 147, 151, 154, We calculated µ C using 1.7 million and 3.9 million heterozygous sites 
159, 160, 165, 179, 182, 186, 187, 189, 200, 211, 219, 225, 228, for the ‘akiapōlā‘au and ‘amakihi, respectively. We used TC values of 
253, 291, 295, 310, 319, 327, 353, 356, 357, 362, 363, 397, 408, 0.43 Mya (Carson and Clague 1995; Fleischer et al. 1998) and 1.0 
409, 412, 422, 432, 453, 466, 468, 491, 493, 513, 541, 634, 638, Mya (Lipman and Calvert 2013).
1020, 1150) scaffolds. The first model may overestimate the long-term genomic substi-
For each of the 3 data sets (the de novo ‘akiapōlā‘au-aligned data tution rate due to the effects of time dependency on molecular rate 
set, the Callicrate-‘amakihi-aligned data set, and the zebra-finch- estimates at recent time scales (e.g., lineage extinction and loss of 
aligned data set), we then calculated autosomal heterozygosities in nonfixed alleles: Ho et al. 2005, 2007) and by ignoring incomplete 
10 000 nonoverlapping windows using VCFtools 0.1.15 (Danecek lineage sorting. To estimate long-term genomic substitution rates, we 
et  al. 2011) following Callicrate et  al. (2014). To identify regions calculated the rate of divergence between the 2 species using the co-
of maximal autosomal divergence between the ‘akiapōlā‘au and alescent formula:
‘amakihi, we calculated densities of fixed autosomal SNVs separ-
ating the 2 species in 10  000  bp nonoverlapping windows using µD = N/L/TMRCA/2
VCFtools 0.1.15 (Danecek et al. 2011). We considered divergence 
windows with a z-score > 7 as candidate adaptive loci for further in- Under this model, µ D is the mean substitution rate based on the di-
vestigation. The stringent z-score cutoff minimized false positives due vergence between the 2 species,  N is the number of single-nucleotide 
to misalignments, sequence and variant errors, and multiple testing substitutions, and TMRCA is the time of most recent common an-
(~95 000 windows per analysis). In the ‘akiapōlā‘au-aligned data set, cestor. We calculated µ D using the total number (9.9 million) of ob-
we discarded contigs and scaffolds below 10 000 bp in length from served substitutions among lineages. The ‘akiapōlā‘au and Hawai‘i 
the heterozygosity and divergence analyses since these had artificially ‘amakihi divergence time was 3.0 Mya (Lerner et  al. 2011). This 
reduced values. Additionally, the ‘amakihi and zebra finch reference model calculates the average substitution rate and does not account 
sequences include artificial “chromosomes” (‘amakihi: Un and Un2; for rate heterogeneity between the honeycreeper branches.
zebra finch: Un) consisting of randomly ordered, concatenated un- Finally, we repeated the divergence-based analysis using the 
placed contigs. We also excluded these “chromosomes” from the het- alignment against the zebra finch. We used 75 million observed sub-
erozygosity and divergence analyses since they had aberrant values stitutions from the zebra finch sequence and an estimated genome 
reflecting misalignments and genotyping errors (data not shown). divergence time from zebra finch of 15 Mya (based on the estimated 
Finally, for each data set, we searched for runs of homozygosity of at separation of Fringillidae [represented by Fringilla] and Estrildidae 
least 1 Mbp using BCFtools 1.9 roh (Li et al. 2009) assuming allele [represented by Estrilda]: Moyle et al. 2016). As this model calcu-
frequencies of 50% (option --Af-dflt 0.50). lates the average rate between the Fringillidae and Estrildidae, the 
estimated honeycreeper rate will only be valid if evolutionary rates 
Estimation of Genomic Substitution Rates are similar between these 2 deeply diverged lineages.
Accurate estimation of genomic substitution rates is challenging 
due to a wide range of factors, including (but not limited to) ac- Genomic Population History Analysis
curacy and precision of selected calibration points (e.g., Fleischer We reconstructed genomic population histories for the 2 species using 
et al. 1998), genomic variant sequencing depth and quality, and time PSMC (Li and Durbin 2011). PSMC estimates ancestral changes in 
dependency of calculated rates (Ho et al. 2005, 2007). Therefore, Ne from a single genome using the coalescent (Li and Durbin 2011). 
we calculated autosomal substitution rates using 3 procedures to The program divides the genome into windows and scores each 
observe the range of estimates. To account for differing levels of window as heterozygous or homozygous. Loci are defined as con-
calibration precision between estimators, we rounded all values to tiguous segments of windows with a similar heterozygosity rate. 
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
Journal of Heredity, 2020, Vol. 111, No. 1 25
PSMC uses a Hidden Markov Model to infer recombination points ‘amakihi identified in the ‘akiapōlā‘au-aligned divergence analysis. 
between loci where the heterozygosity rate changes. The number For comparison, we investigated Ensembl-annotated zebra finch 
of heterozygous windows within each locus is proportional to the genes (annotation taeGut3.2.4) that overlapped divergent regions in 
TMRCA of that window. The number of loci at each TMRCA is inversely the ‘amakihi-aligned and zebra-finch-aligned data sets.
proportional to Ne at that time. The resulting population history is 
then scaled by the mutation rate and the generation time. The in- Bill Morphology Gene Analysis
terpretation of PSMC is thus sensitive to inaccurate mutation rates We compiled a list of genes previously demonstrated to be as-
and generation times (Li and Durbin 2011). PSMC does not recover sociated with bill morphology in birds (ALX1, BMP4, CALM1, 
recent events accurately due to the limited number of heterozygous CALM2, CALML3, DLK1, HMGA2: Abzhanov et al. 2004, 2006; 
SNVs. It also tends to interpret rapid population size changes (e.g., Lamichhaney et  al. 2015, 2016; Chaves et  al. 2016). We aligned 
a strong bottleneck effect) as gradual clines and can be confused by these candidate genes’ protein sequences annotated in Ensembl for 
population structure (Li and Durbin 2011). zebra finch (assembly taeGut 3.2.4; Warren et al. 2010) against the 
To control for the impact of reference bias on demographic re- ‘akiapōlā‘au and ‘amakihi genome sequences using TBLASTN 2.7.0+ 
construction, we analyzed the de novo ‘akiapōlā‘au genome aligned (maximum expect value: 0.1; Altschul et al. 1997). We retained the 
sequences, the Callicrate ‘amakihi assembly aligned sequences, and best hits for each exon. We then manually inspected retained align-
the zebra-finch-aligned genomes independently. We explored the im- ments to identify sequence variants for further analysis. To control 
pact of the concatenated, unplaced contig “chromosomes” (Un for for reference bias in the ‘amakihi genome, we assessed both the 
zebra finch; Un and Un2 for ‘amakihi) on demographic reconstruc- Callicrate assembly and the ‘akiapōlā‘au-realigned sequence.
tions by repeating the analyses with these sequences excluded. As TBLASTN did not identify exon 1 of HMGA2 in both the 
recommended by the PSMC authors, we generated depth-filtered ‘akiapōlā‘au and ‘amakihi genomes (see below). To confirm this re-
consensus sequences for each assembly setting the minimum and sult, we searched for the corresponding exon 1 DNA sequence using 
maximum depths to one-third and twice the average short-insert megaBLAST (BLASTN 2.7.1+: Camacho et al. 2008).
sequencing depths (38×–228× for ‘akiapōlā‘au and 20×–120× for Since we identified 2 DLK1 amino acids that differed between 
‘amakihi) using BCFtools 1.9 mpileup and BCFtools call (option -c) the ‘akiapōlā‘au and ‘amakihi (see below), we investigated this gene 
(Li et al. 2009). Variant calls were converted to FASTQ format using phylogenetically to determine the uniqueness of these variants. We 
vcfutils.pl vcf2fq. We retained variants with qualities of at least 20 obtained DLK1 protein sequences predicted using the National 
for the analysis. Bootstrapped PSMC (100 replicates) was performed Center for Biotechnology Information annotation pipeline for 66 
using PSMC 0.6.5 (options -N 25, -t 15, -r 5, -p “64*1”). We scaled avian species in GenBank. We aligned them using the Geneious 
the PSMC reconstructions using the range of estimated genomic sub- and ClustalW 2.1 (Larkin et al. 2007) aligners in Geneious Prime 
stitution rates (see above) and generation times of 2 years for the 2019.0.4 (Biomatters Ltd, Auckland, New Zealand). We trimmed 
‘akiapōlā‘au and 1 year for the ‘amakihi (Pratt 2005, p. 154). We the alignment to the Ensembl-annotated zebra finch sequence. We 
also scaled the PSMC reconstruction using the collared flycatcher manually inspected the trimmed alignment for erroneous stop co-
(Ficedula albicollis) germ-line mutation rate (4.6  × 10−9 substitu- dons and indel misalignments. We then constructed a rapid boot-
tions/site/generation; Smeds et  al. 2016) since no drepanid rate is strapping (100 replicates), maximum likelihood phylogenetic tree 
currently available. using the RAxML 8.2.11 under the GAMMA BLOSUM62 protein 
substitution model with a random starting tree (Stamatakis 2014).
Genome Annotation and Gene Prediction
We annotated the de novo ‘akiapōlā‘au assembly, the realigned Results
‘amakihi consensus sequence, and the original Callicrate ‘amakihi 
genome sequence for comparison. We trained AUGUSTUS gene pre- Genome Assembly
diction using the BUSCO 3.0.2 retraining results (Waterhouse et al. We sequenced the ‘akiapōlā‘au genome to a nominal coverage of 
2017). We annotated repeats for the 3 assemblies using RepeatMasker ~134× (Table 1). We generated 915 million short-insert paired reads 
4.0.7 (Smit et al. 2013–2015) and Repbase Update 20170127 (all (~114×), 107 million 3–5 kbp mate pairs (~13×), and 58 million 8–10 
organisms; Bao et al 2015). We obtained the zebra finch (assembly kbp (~7×). The genome was assembled into 40 124 contigs (contig 
taeGut 3.2.4; Warren et  al. 2010) proteins and cDNAs from the N50: 75 423 bp) and 9908 scaffolds (scaffold N50: 3 340 285 bp). 
Ensembl database. We mapped the zebra finch proteins and cDNAs A  total of 32 612 (98.7%) contigs were assembled into scaffolds. 
to the genome assemblies using Exonerate 2.2 protein2genome The total contig assembly length was 1  014  901  717  bp and the 
(Slater and Birney 2005) and BLAT 36x1 (minimum sequence iden- total scaffolded length was 1 032 508 528 bp, similar to previously 
tity 92%; Kent 2002), respectively. We converted the annotated sequenced bird genomes (Jarvis et al. 2014).
repeats, proteins, and cDNAs to hints for AUGUSTUS 3.3 and pre- Genic content was very complete for all assemblies. BUSCO 
dicted genes using the BUSCO-trained model (Stanke et al. 2006a, found ~94% of the 4915 BUSCO groups in all assemblies, with low 
2006b, 2008). Hint weights used AUGUSTUS defaults except that rates of duplication (~1%) and missingness (2% in the ‘akiapōlā‘au, 
the RepeatMasker hints were weighted as 1.15 for “nonexonpart” 2.5% in the ‘amakihi) (Table 2). We observed slightly better BUSCO 
and 1 for all other states. We identified orthologous predicted genes scores for the original Callicrate ‘amakihi assembly than the re-
between the 3 assemblies using Proteinortho 5.16b (Lechner et al. aligned sequence, reflecting variance between BUSCO runs or pos-
2011). We classified predicted proteins’ domains, functions, and sible reference bias.
pathways using InterProScan 5.33-72.0 (options -goterms -pa; 
Jones et  al. 2014) against all eukaryotic databases. We used the Genomic Diversity and Divergence
InterProScan classifications to investigate predicted proteins that In the ‘akiapōlā‘au-aligned analyses, we observed 1 665 867 hetero-
overlapped highly divergent regions between the ‘akiapōlā‘au and zygous autosomal SNVs in the ‘akiapōlā‘au or one every 570 bases 
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
26 Journal of Heredity, 2020, Vol. 111, No. 1
Table 1. ‘Akiapōlā‘au genome assembly statistics
Total scaffold length: 1 032 508 528 bp Total contig length: 1 014 901 717 bp
Number scaffolds: 9908 Number contigs: 40 124
Scaffold N50: 3 340 285 bp Contig N50: 75 423 bp
Scaffold L50: 85 Contig L50: 3815
Longest scaffold: 15 969 922 bp Longest contig: 629 073 bp
Mean scaffold length: 104 210 bp Mean contig length: 25 294 bp
Median scaffold length:: 1321 bp Median contig length: 7069 bp
Contig %A: 29.03% Contig %C: 20.95%
Contig %G: 20.97% Contig %T: 29.04%
Mean break length between scaffolded contigs 582 bp Nominal sequencing coverage: 134×
Table 2. BUSCO results local heterozygosity, especially at the telomeres of the chromo-
somes, suggestive of misaligned, collapsed repetitive elements 
BUSCO  ‘Akiapōlā‘au ‘Amakihi ‘Amakihi  
(Supplementary Figures 1 and 2; Treangen and Salzberg 2012). 
category (Callicrate) (realigned)
This indicates that these 2 data sets include less accurately assem-
Complete 4635 (94.3%) 4609 (93.8%) 4604 (93.7%) bled regions excluded from the de novo ‘akiapōlā‘au genome as-
Complete,  4575 (93.1%) 4560 (92.8%) 4549 (92.6%) sembly, decreasing overall analytical quality.
single-copy While local genomic divergences were relatively constant across 
Complete,  60 (1.2%) 49 (1.0%) 55 (1.1%) the genome within all 3 data sets (‘akiapōlā‘au-aligned: 0.00439 ± 
duplicated 0.00143 fixed SNVs/bp; ‘amakihi-aligned: 0.00429 ± 0.00170 fixed 
Fragmented 181 (3.7%) 181 (3.7%) 189 (3.8%) SNVs/bp; zebra-finch-aligned: 0.00321  ± 0.00118 fixed SNVs/bp), 
Missing 99 (2.0%) 125 (2.5%) 122 (2.5%) there was decreased mean divergence with reference distance from 
Total BUSCOs 4915 4915 4915
‘akiapōlā‘au. In the ‘akiapōlā‘au-aligned data set, divergence ex-
searched
ceeded 0.015 fixed SNVs/bp (z-scores > 7.4) in 7 distinct regions on 
‘akiapōlā‘au scaffolds 39, 62, 119, 349, 548, 788, and 810 (Figure 
Total numbers of identified BUSCOs of each category are given, with the 
percentage of the total found in parentheses. 2). The elevated diversity regions on scaffolds 39, 62, and 119 were 
50–70 kbp long, while the remaining were 10–40 kbp. We observed 
across the ~950 Mbp autosomal genome. The realigned ‘amakihi had 20 divergent regions in the ‘amakihi-aligned data set across chromo-
3 938 436 or one every 241 bases. The 2 genomes were separated somes 1 (n = 1), 2 (n = 5), 3 (n = 1), 4 (n = 2), 4A (n = 1), 8_random 
by 4 174 597 fixed SNVs. The ‘akiapōlā‘au’s local heterozygosities (n  = 1), 9 (n  = 1), 10 (n  = 1), 15 (n  = 2), 18 (n  = 1), 20 (n  = 1), 
(mean ± standard deviation: 0.00172 ± 0.00111 heterozygous sites/ 26 (n = 1), 22_random (n = 1), and 27 (n = 1). Most (14) of these 
bp) were less than half (0.416×; unpaired t-test P << 0.0001) the were ~10 kbp long. Six longer (30–60 kbp) divergent regions were 
‘amakihi’s (0.00413 ± 0.00192 heterozygous sites/bp) (Figure 2). We observed on chromosomes 4, 4A, 8_random, 9, 10, and 15. We ob-
observed one 590 kbp region of the ‘amakihi genome (corresponding served only 4 divergent regions in the zebra-finch-aligned data set: 3 
to ‘akiapōlā‘au scaffold 261) with exceptionally high heterozygosity of these matched the longer divergent regions on chromosomes 4A, 
(0.0127 heterozygous sites/bp, z-score = 4.46), suggestive of a col- 8_random, and 10 in the ‘amakihi-aligned data set and the remaining 
lapsed repetitive region unique to the ‘amakihi genome (Treangen was a short (~10 kbp) region on chromosome 2 close (~55 kbp away) 
and Salzberg 2012). We also observed elevated heterozygosities in to one of the divergent regions in the ‘amakihi-aligned data set.
the scaffolds between 10 and 100 kbp (especially in the ‘akiapōlā‘au), We detected no long (1 Mbp) runs of homozygosity in either 
which likely derive from assembly or alignment errors of repetitive genome in the ‘akiapōlā‘au- or zebra-finch-aligned sequences. In 
sequences (Treangen and Salzberg 2012). agreement with Callicrate et  al. (2014), we observed 8 long runs 
There were 1 705 302 ‘akiapōlā‘au SNVs and 3 856 230 ‘amakihi of homozygosity on chromosomes 1 (n  = 2) and 6 (n  = 6) in the 
heterozygous autosomal SNVs in the ‘amakihi-aligned data set ‘amakihi in the ‘amakihi-aligned data set. No runs of homozygosity 
and 1  931  190  ‘akiapōlā‘au SNVs and 3  063  769  ‘amakihi het- were detected in the ‘akiapōlā‘au in this data set. The absence of 
erozygous autosomal SNVs in the zebra-finch-aligned data set. these runs of homozygosity in the other data sets are likely due to 
In the ‘amakihi-aligned data set, local heterozygosities were a combination of inaccurate allele calls in the divergent reference 
similar to the ‘akiapōlā‘au-aligned data set in the ‘akiapōlā‘au sequences and the absence of chromosome-level contiguity informa-
(0.00174 ± 0.00159 heterozygous sites/bp), but slightly decreased tion in the scaffold-level ‘akiapōlā‘au genome.
in the ‘amakihi (0.00393  ± 0.00212 heterozygous sites/bp). The Overall, we found strong evidence of reference bias impacting 
zebra-finch-aligned data set inferred higher ‘akiapōlā‘au local estimates of heterozygosity and divergence and confounding iden-
heterozygosities (0.00197  ± 0.00194 heterozygous sites/bp) and tification of runs of homozygosity (Supplementary Figures 1 and 
substantially lower ‘amakihi values (0.00312  ± 0.00198 hetero- 2). As expected, alignment against the ‘amakihi genome produced 
zygous sites/bp). As a result the ratio of ‘akiapōlā‘au to ‘amakihi heterozygosity and divergence estimates that were more similar to 
heterozygosities increased in these data sets (‘amakihi-aligned: those generated using the ‘akiapōlā‘au as reference than either were 
0.443×; zebra-finch-aligned: 0.631×; unpaired t-test P << 0.0001 in to those using the more divergent zebra finch.
both cases). Variance in local heterozygosities was greater in both 
the ‘amakihi-aligned and zebra-finch-aligned data sets than in the Genomic Substitution Rates
‘akiapōlā‘au-aligned one. In the ‘amakihi-aligned and zebra-finch- Under the first model (µ C), the ‘akiapōlā‘au rate was 0.9–2.1 × 10
−9 sub-
aligned data sets, both species had numerous peaks of apparent stitutions/site/year, and the ‘amakihi µ C was 2.1–4.8 × 10
−9 substitutions/
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
Journal of Heredity, 2020, Vol. 111, No. 1 27
Figure 2. Genomic heterozygosities and SNV densities of the ‘akiapōlā‘au and ‘amakihi calculated in 100 kbp windows. Colors correspond to individual scaffolds. 
Scaffolds are sorted by decreasing length with longer ones to the left. The de novo ‘akiapōlā‘au genome was used as the reference sequence.
site/year. The increased rate in the ‘amakihi likely reflects its higher ef- finch population is large (IUCN 2019), permitting the evolution and 
fective population size and/or decreased rates of purifying selection due maintenance of an increased number of mutations.
to being an ecological generalist. We estimated the mean divergence 
substitution rate (µ D) as 1.7 × 10
−9 substitutions/site/year (divergence Genomic Population History Analysis
between ‘akiapōlā‘au and ‘amakihi) to 2.6 × 10−9 substitutions/site/year The PSMC reconstructions showed markedly divergent population 
(divergence from zebra finch). The increased substitution rate from the histories for the 2 species (Figure 3). The ‘amakihi shows an ini-
zebra finch-based estimate compared to the within-honeycreeper esti- tial decrease in Ne ending 2–5 Mya in the ‘akiapōlā‘au-aligned and 
mates could represent a faster evolutionary rate along the Estrildid lin- 1–3 Mya in the zebra-finch-aligned data sets, respectively (Figure 3; 
eage artificially inflating the estimated honeycreeper rate. An increased Supplementary Figures 3 and 4). We did not detect this decline in 
evolutionary rate along the Estrildid lineage is likely since the zebra the ‘amakihi aligned data set. After a period of stability, the ‘amakihi 
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
28 Journal of Heredity, 2020, Vol. 111, No. 1
500 rate scaling (Supplementary Figure 4). Exclusion of the artificial 
‘Akiapōlā‘au
‘Amakihi “chromosomes” Un and Un2 had negligible impact on demographic 
400 reconstructions using the ‘amakihi reference (Supplementary Figure 
5). In the zebra-finch-aligned data set, exclusion of “chromosome” 
Un had minimal impact on the ‘akiapōlā‘au reconstruction, but re-
300
duced inferred Ne in the ‘amakihi (Supplementary Figure 6). The 
‘amakihi’s demographic history topology however remained similar.
200
Genome Annotation and Candidate Gene 
100
Identification
AUGUSTUS predicted 34  450 genes in the ‘akiapōlā‘au genome, 
0
104 105 106 107 108 33 083 genes in the Callicrate ‘amakihi assembly, and 34 620 genes 
Years before present in the realigned ‘amakihi sequence. Proteinortho identified 27 774 
orthologous clusters between the 3 assemblies. Thirteen predicted 
Figure 3. Pairwise sequential Markovian coalescent reconstruction of ‘akiapōlā‘au and 10 predicted ‘amakihi genes overlapped the highly 
‘akiapōlā‘au and ‘amakihi demographic histories. Demographic reconstructions divergent genomic regions. One gene classified as “hydrocephalus-
used substitution rates assuming colonization of Hawai‘i island 1.0 Mya. inducing-like” with immunoglobulin-like folds was found in both 
The ‘akiapōlā‘au and ‘amakihi generation times were 2  years and 1  year, genome sequences on ‘akiapōlā‘au scaffold 788. A  cadherin-like 
respectively. Initial results are plotted using darker lines, with bootstrap 
gene (scaffold 119) and an ACE1-Sec16-like gene (scaffold 62) were 
replicates in lighter hues.
identified in the ‘akiapōlā‘au. We found 5 genes deriving from retro-
lineage underwent a rapid growth phase starting 0.3–1 Mya in the viruses: 3 in the ‘amakihi on scaffold 119 (an avian retrovirus enve-
‘akiapōlā‘au-aligned and ‘amakihi-aligned data sets and 0.2–0.8 lope protein, an integrase, and an RNase H), another RNase H in the 
Mya in the zebra-finch-aligned data set. In the ‘akiapōlā‘au-aligned ‘amakihi on scaffold 548, and a retroviral aspartyl protease in the 
and ‘amakihi-aligned data sets, this was followed by a secondary ‘akiapōlā‘au on scaffold 548. The remaining predicted genes had no 
population decline starting 70–200 thousand years ago (kya) and identified functions.
ending between 40–90 kya (‘akiapōlā‘au-aligned data set) or 20–70 The divergent regions from the ‘amakihi-aligned data set 
kya (‘amakihi-aligned data set). In both reconstructions, this was overlapped 10 protein coding genes (ENSTGUG00000000646, 
followed by a period of population stability or slight expansion until ENSTGUG00000001862, ENSTGUG00000002283, ENSTGUG 
10–30 kya. Conversely, in the zebra-finch-aligned reconstructions, 00000009266, FOXB1, KIF6, MAML3, PHF20, RBP1, and 
instead of a decline after the expansion, there is a ~50 thousand-year TIMM17A) and an RNA gene (ENSTGUG00000018063). Only 
period of population stability. The stable period is followed by a sec- ENSTGUG00000009266, FOXB1, and ENSTGUG00000018063 
ondary recent expansion (Supplementary Figure 3). Depending on were found using the zebra-finch-aligned data set.
the substitution rate, we detected a secondary rapid recent decline 
in the ‘akiapōlā‘au-aligned and the zebra-finch-aligned ‘amakihi data Bill Morphology Genes
sets. As the confidence interval around this secondary decline is quite All candidate bill morphology genes were discovered on single 
large, this is likely an artifact due to there being insufficient hetero- ‘akiapōlā‘au scaffolds, supporting the accuracy of the ‘akiapōlā‘au 
zygous SNVs to infer more recent demographic history accurately. genome assembly. These genes were also identified in the original 
The ‘akiapōlā‘au N  remained relatively stable, but smaller than Callicrate ‘amakihi assembly to be located on the expected chromo-e
the ‘amakihi’s, throughout the population history since after the somes based on assumed conserved synteny with the zebra finch. 
‘amakihi population expansion. The ‘akiapōlā‘au underwent slow This indicates that we recovered the correct orthologs for these 
population growth after ~1 Mya, with a recent population decline proteins.
starting 15–50 kya. Additionally, using the substitution rate derived We recovered the coding sequence of ALX1 (4 exons), BMP4 (2 
from the emergence of Hawai‘i at 1 Mya and the zebra finch ref- exons), and CALML3 (1 exon) in their entirety. The 3 genes were lo-
erence, the ‘akiapōlā‘au undergoes an apparent dramatic expansion cated on ‘akiapōlā‘au scaffolds 83, 436, and 751 (Callicrate chromo-
within the last ~20 000 years (Supplementary Figure 3). This is likely somes 1A, 5, and 1A), respectively. The ‘akiapōlā‘au and ‘amakihi 
an artifact due to insufficient variable sites at recent time scales. The amino acid sequences were identical for all 3 genes. The honey-
confidence interval includes population stasis to expansion. The ap- creeper ALX1 differed from the zebra finch peptide by 3 amino acid 
pearance of this expansion is dependent on the substitution rate and replacements (p.L103F, p.I221V, and p.T266S). The honeycreeper 
does not appear with higher rates (data not shown). BMP4 differed from the zebra finch’s by 2 amino acid replacements 
While the reconstructions using the ‘akiapōlā‘au and ‘amakihi (p.H23Y and p.T213S). CALML3 differed by one replacement 
references were similar to each other, we observed significant im- (p.K31R).
pact on inferred demographic history using the zebra finch refer- CALM1 and CALM2 (both calmodulins) produce identical 
ence. ‘Amakihi Ne estimates and degrees of population size changes amino acid sequences. We identified both copies of this sequence on 
were markedly smaller in the zebra-finch-aligned reconstructions ‘akiapōlā‘au scaffolds 54 and 192. The scaffold 54 sequence corres-
(Supplementary Figure 3). The decreased ‘amakihi population sizes ponded to Callicrate ‘amakihi chromosome 3, while the scaffold 192 
and apparent population stasis period are likely derived from in- sequence matched Callicrate ‘amakihi chromosome 5. This indicates 
correct variant calls in the zebra-finch-aligned reconstructions since that the scaffold 54 copy is CALM2 and the scaffold 192 sequence is 
they are not concordant with the ‘akiapōlā‘au-aligned or ‘amakihi- CALM1 assuming that the original ‘amakihi assembly was correctly 
aligned results. Moreover, N  estimates for both species and timing oriented to the zebra finch in this region. Nevertheless, both copies e
of the population size changes were determined by the substitution covered amino acids 11 to 139 of the calmodulin sequences (from an 
Effective population size (× 104)
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
Journal of Heredity, 2020, Vol. 111, No. 1 29
expected length of 148 amino acids across 6 exons). No amino acid before its ancestors’ colonization of Hawai‘i island. We posit that 
replacements were observed in either the ‘akiapōlā‘au or ‘amakihi. that the ‘amakihi rapid population growth period corresponds to 
While the initial TBLASTN search only identified exon 2 of the colonization of Hawai‘i and subsequent expansion into new 
HMGA2 (from a total of 2 exons) on ‘akiapōlā‘au scaffold 24 habitat. Scaling the reconstruction using most estimated rates dates 
(Callicrate ‘amakihi chromosome 1A), the subsequent megaBLAST the onset of the expansion to 0.5–1 Mya (Figure 3; Supplementary 
search found exon 1 on the same scaffold as exon 2. The honey- Figure 3), correlating well with more recent estimates of the first sub-
creeper peptide sequences were identical and shared a single replace- aerial emergence of the island (Lipman and Calvert 2013). The latest 
ment (p.A33T) from the zebra finch peptide. estimated dates for the start of the growth phase is ~200–300 kya 
For DLK1, we recovered all 4 exons in their entirety for the (depending on reference sequence). However, these dates derive from 
‘akiapōlā‘au and the realigned ‘amakihi (‘akiapōlā‘au scaffold 27/ the substitution rates estimated with the emergence of the youngest 
Callicrate ‘amakihi chromosome 5). Exon 1 was missing from the island of Hawai‘i at 0.43 Mya, which is probably an underestimate 
original ‘amakihi genome assembly. Phylogenetic analysis clustered of its age, and the flycatcher germ-line mutation rate, which may 
the honeycreeper DLK1 sequences with canary (Serinus canaria: be inappropriate for honeycreepers. Nevertheless, in all cases, the 
GenBank accession XP_009099991.2; Supplementary Figure 7; del ‘amakihi expansion followed shortly after the date used to anchor 
Hoyo et al. 2019). The ‘akiapōlā‘au, ‘amakihi and canary shared 5 the emergence of Hawai‘i, strongly supporting the argument that 
amino acid replacements from the zebra finch (p.N3S, p.Q106K, this expansion represents a colonization event. The slow population 
p.Y176H, p.T188P, and p.K223R). The ‘amakihi had an asparagine growth in the ‘akiapōlā‘au after ~1 Mya may also represent its col-
at amino acid 76, while the ‘akiapōlā‘au, canary, and zebra finch onization of Hawai‘ and/or population expansion due to bill mor-
shared a threonine residue at this position. The ‘amakihi asparagine phological adaptation permitting better exploitation of its dietary 
residue was unique across all investigated avian species. At amino niche compared to a presumed nukupu‘u-like ancestor.
acid 220, the ‘akiapōlā‘au and canary had a valine residue, while the The recent population declines in both species are more difficult 
zebra finch had a methionine. Amino acid 220 was not resolved in to explain. One possibility is that it may represent the recent 19th 
the realigned ‘amakihi, but was an alanine in the original assembly. and 20th century bottleneck in the honeycreeper populations due to 
Inspection of the ‘amakihi realignment revealed a C/T SNV at amino the introduction of Avipoxvirus and avian malaria to the Hawai‘ian 
acid 220, indicating that both alanine and valine residues were pre- islands (Warner 1968; van Riper et  al. 1986; Samuel et  al. 2011, 
sent at this site. 2015). The sequenced ‘amakihi was a low elevation bird (Callicrate 
et  al. 2014; Cassin-Sackett et  al. 2019), so the bottleneck due to 
the introduction of avian malaria and the strong selection for toler-
Discussion ance to the parasite is probably reflected in its genome (Foster et al. 
Genomic Diversity and Demography 2007; Atkinson et al. 2013). PSMC often incorrectly infers gradual 
Our heterozygosity estimates and demographic reconstructions sup- declines instead of bottlenecks and struggles to accurately infer re-
ported the expected genomic patterns for the 2 Hawai‘ian honey- cent demographic events due to limited variation and recombination 
creeper species. The ‘akiapōlā‘au is less heterozygous and has a events (Li and Durbin 2011). Alternatively, these may represent true 
smaller effective population size than the ‘amakihi, consistent with population declines in the past. For instance, population declines 
an ecological specialist restricted to a small endemic range with low could derive from restriction of suitable habitat caused by climate 
census population numbers. The ‘amakihi is more heterozygous change (e.g., increased global temperatures since the Last Glacial 
and has a much greater N , consistent with its comparatively larger Maximum). Further research using specimens that predate the re-e
census population and wider ecological niche. As previously reported cent bottleneck will clarify the timing of this population decline. 
by Callicrate et al. (2014), the ‘amakihi’s runs of homozygosity on Similarly, analysis of Hawai‘i ‘amakihi individuals from populations 
chromosomes 1 and 6 are suggestive of a selective sweep on the same presumably less affected by malaria and the bottleneck (e.g., from 
chromosomes, possibly due to the recent evolution of tolerance for high elevations) could help determine the population specificity and 
avian malaria in lowland ‘amakihi. Encouragingly for its long-term timing of the observed decline.
survival, the ‘akiapōlā‘au’s individual heterozygosity (~0.0017 het-
erozygous sites/bp) is higher than the range previously observed for Honeycreeper Adaptive Radiation Genomics
threatened bird species (0.0001–0.00091 heterozygous sites/bp: Li Our genomic scan identified small strongly divergent regions that 
et al. 2014; Cortes-Rodriguez et al. 2019) such as New Zealand hihi distinguished the ‘akiapōlā‘au and ‘amakihi genomes, suggesting that 
(Notiomystis cincta; 0.00069 heterozygous sites/bp: de Villemereuil the honeycreeper radiation is driven by diversifying selection on a 
et  al. 2019), Hawai‘ian ‘alala (Corvus hawaiiensis; 0.0004 sites/ small number of genes and pathways rather than uniform selection 
bp: Sutton et al. 2018), and Marianas åga (Corvus kubaryi; 0.0001 across the genome (Chaves et  al. 2016; Wolf and Ellegren 2017). 
sites/bp: Cortes-Rodriguez et  al. 2019). Furthermore, the absence This pattern has been observed in birds previously: for instance, only 
of large runs of homozygosity in the ‘akiapōlā‘au suggests that it is 6 genomic regions differed between golden-winged and blue-winged 
not recently inbred (e.g., Dobrynin et al. 2016). Additional genome warblers, of which 4 were linked to plumage genes (Toews et  al. 
sequences are needed to determine whether inbreeding is occurring 2016). While the honeycreeper genomic islands of differentiation 
in other currently unsampled individuals. may contain adaptation or speciation genes, further analysis is re-
Despite the current range of substitution rate estimates, our quired to determine whether the identified genes in these regions are 
PSMC reconstructions permit us to tentatively identify genomic truly involved in honeycreeper speciation or are the result of other 
events likely related to the honeycreeper radiation. While we are un- processes (e.g., linked selection or population divergence; Wolf and 
able to link the initial population decline in the ‘amakihi to a known Ellegren 2017). Nevertheless, the identified genes in these regions are 
historical event due to its imprecise dating (ending 1–5 Mya), it may candidates to explain honeycreeper evolution. Intriguingly, FOXB1 
represent a bottleneck caused by a colonization or speciation event expression regulates numerous developmental processes, notably 
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
30 Journal of Heredity, 2020, Vol. 111, No. 1
neurological development and visual learning (UniProt Consortium evolution within DLK1, it does not rule out variation in DLK1 ex-
2019). Evolution in FOXB1 expression could contribute to behavior pression as playing a role.
differences between the 2 species. Furthermore, in humans, genic An obvious next step to better understand the evolution of 
regulation of brain growth and development is involved with the the ‘akiapōlā‘au’s unique bill morphology would be to search the 
development of congenital hydrocephalus (a disorder characterized nukupu‘us’ genomes for signatures of selection that could explain 
by abnormal accumulation of cerebral spinal fluid within the cere- this phenomenon. Due to the nukupu‘us’ apparent extinction, such 
bral ventricles; Kahle et  al. 2016). Therefore, similar to FOXB1, an analysis will require the use of museum specimens and ancient 
the identified “hydrocephalus-inducing-like” gene suggests genic DNA techniques. In addition, developmental transcriptomic ana-
regulation of neurological development in honeycreeper evolution. lyses of Hawaiian honeycreepers that are not endangered with 
The predicted ‘akiapōlā‘au cadherin-like gene suggests a role in differing bill morphologies (as has been done for Darwin’s finches: 
the calcium-binding pathway in honeycreeper evolution. The cal- Abzhanov et al. 2004, 2006) may enable discovery of other candi-
modulin pathway (which also involves calcium binding) was previ- date genes involved in the differentiation of bill morphology in this 
ously shown to influence bill length in Darwin’s finches (Abzhanov extensive adaptive radiation.
et  al. 2006). It is therefore plausible that this gene is involved in 
the ‘akiapōlā‘au’s unique bill morphology. While the identified 
avian retroviruses may also play a role in honeycreeper evolution, Conclusions
endogenous retrovirus expression and regulation remains poorly Our de novo ‘akiapōlā‘au genome assembly provides a critical resource 
understood (Hu et al. 2016). for genomic investigation of the endangered honeycreeper radiation. 
The remaining named identified genes are involved in a variety Comparison of the ‘akiapōlā‘au and ‘amakihi genomes revealed the ex-
of cellular processes: KIF6 is involved in ATP binding and micro- pected pattern of reduced heterozygosity and effective population size 
tubule movement, MAML3 is part of the Notch signaling pathway in the ‘akiapōlā‘au. The age of the island on which both species occur 
and regulates transcription, PHF20 contributes to histone acetyl- provides, with assumptions, a means to assess comparative substitution 
ation, DNA binding, and transcription regulation, RBP1 binds ret- rates within each species. Speciation within the drepanids appears to be 
inol, and TIMM17A regulates transmembrane protein transport driven by selection on a small number of islands of divergence, rather 
(UniProt Consortium 2019). The ‘akiapōlā‘au ACE1-Sec16-like than uniform diversification across the genome. Most notably, the gene 
gene is too poorly identified to reliably infer function, but is likely in- FOXB1 was discovered within one of these genomic islands of diver-
volved in protein transport or binding (UniProt Consortium 2019). gence, suggesting a role in neurological or behavioral adaptation be-
Due to the low-level nature of these genes, selection on them likely tween the 2 species. The amino acid changes in only one bill morphology 
contributes to a range of phenotypes, which cannot be disentangled associate gene, DLK1, indicates that it could be a strong candidate for 
without further controlled experimental study. The zebra-finch- controlling bill morphology differences in honeycreepers.
predicted (ENSTGUG00000000646, ENSTGUG00000001862, 
ENSTGUG00000018063, ENSTGUG00000002283, 
ENSTGUG00000009266) and remaining ‘akiapōlā‘au/‘amakihi- Supplementary Material
predicted genes’ functions are currently unknown, but could be 
Supplementary data are available at Journal of Heredity online.
critical to honeycreeper evolution.
With the exception of DLK1, we found no evidence of between-
honeycreeper differences in the candidate bill morphology genes de- Data Availability
rived from the literature, suggesting adaptation in different genes or 
in regulatory regions (Wittkopp and Kalay 2012). Gene regulation, The ‘akiapōlā‘au genome sequence and raw data are available 
rather than protein sequence evolution, can be a primary source of in GenBank (accession SGIP00000000) and as NCBI BioSample 
phenotypic variation leading to speciation (Mack and Nachman SAMN10867508 under BioProject PRJNA520783.
2017). Gene regulation evolution could explain the concurrent 
rapid phenotypic and limited genotypic evolution within the honey- Funding
creeper clade (Lerner et  al. 2011). Since we observed variation in 
DLK1 coding sequence between honeycreeper species, the gene is a This work was supported by the Smithsonian Institution; and the 
candidate for further analysis to explain honeycreeper bill morph- National Science Foundation (grant numbers 1547168, 1717498).
ology. DLK1 plays roles in calcium ion binding, cell differentiation, 
and negative regulation of the Notch signaling pathway (UniProt 
Consortium 2019). DLK1 has previously been associated with Acknowledgments
variation in bill size in Darwin’s finches (Lamichhaney et al. 2015; We thank Jack Jeffrey for providing the honeycreeper photographs. This re-
Chaves et al. 2016). If honeycreeper bill morphologies evolved using search used the Smithsonian Institution High Performance Cluster (SI/HPC).
the same gene, this would suggest that bill morphology evolution 
is constrained to a limited set of genic pathways within passerines. References
Interestingly, we observed the unique DLK1 mutation in the gener-
alist ‘amakihi, rather than the specialist ‘akiapōlā‘au. This is contrary Abzhanov A, Kuo WP, Hartmann C, Grant BR, Grant PR, Tabin CJ. 2006. 
The calmodulin pathway and evolution of elongated beak morphology in 
to the expectation that the more derived morphology would cor-
Darwin’s finches. Nature. 442:563–567.
respond to the derived gene sequence. However, while not as mor- Abzhanov A, Protas M, Grant BR, Grant PR, Tabin CJ. 2004. Bmp4 and mor-
phologically unique, the ‘amakihi has also evolved a different bill phological variation of beaks in Darwin’s finches. Science. 305:1462–1465.
phenotype from that of its presumed granivorous finch-like ancestor Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, Lipman DJ. 
(Lerner et al. 2011). While this also indicates that the ‘akiapōlā‘au’s 1997. Gapped BLAST and PSI-BLAST: a new generation of protein data-
unique bill morphology does not derive from protein sequence base search programs. Nucleic Acids Res. 25:3389–3402.
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
Journal of Heredity, 2020, Vol. 111, No. 1 31
Atkinson CT, Saili KS, Utzurrum RB, Jarvi SI. 2013. Experimental evidence for chondrial DNA sequence and osteological characters. Stud Avian Biol. 
evolved tolerance to avian malaria in a wild population of low elevation 22:98–103.
Hawai’i ‘amakihi (Hemignathus virens). Ecohealth. 10:366–375. Foster  JT, Woodworth  BL, Eggert  LE, Hart  PJ, Palmer  D, Duffy  DC, 
Bao W, Kojima KK, Kohany O. 2015. Repbase update, a database of repetitive Fleischer RC. 2007. Genetic structure and evolved malaria resistance in 
elements in eukaryotic genomes. Mob DNA. 6:11. Hawaiian honeycreepers. Mol Ecol. 16:4738–4746.
Beadell JS, Ishtiaq F, Covas R, Melo M, Warren BH, Atkinson CT, Bensch S, Futuyma DJ. 1986. Evolutionary biology. Sunderland (MA): Sinauer.
Graves GR, Jhala YV, Peirce MA, et al. 2006. Global phylogeographic Gorresen PM, Camp RJ, Pratt TK. 2007. Forest bird distribution, density and 
limits of Hawaii’s avian malaria. P R Soc B. 273: 2935–2944. trends in the Ka‘ū region of Hawai‘i Island. U.S. Geological Survey Open-
Berner D, Salzburger W. 2015. The genomics of organismal diversification illu- File report 2007-1076. U.S. Geological Survey. Available from: https://
minated by adaptive radiations. Trends Genet. 31:491–499. pubs.usgs.gov/of/2007/1076
Bradnam  KR, Fass  JN, Alexandrov  A, Baranay  P, Bechner  M, Birol  I, Henning F, Meyer A. 2014. The evolutionary genomics of cichlid fishes: ex-
Boisvert S, Chapman JA, Chapuis G, Chikhi R, et al. 2013. Assemblathon plosive speciation and adaptation in the postgenomic era. Annu Rev Gen-
2: evaluating de novo methods of genome assembly in three vertebrate omics Hum Genet. 15:417–441.
species. Gigascience. 2:10. Ho SY, Phillips MJ, Cooper A, Drummond AJ. 2005. Time dependency of mo-
Brawand  D, Wagner  CE, Li  YI, Malinsky  M, Keller  I, Fan  S, Simakov  O, lecular rate estimates and systematic overestimation of recent divergence 
Ng AY, Lim ZW, Bezault E, et al. 2014. The genomic substrate for adap- times. Mol Biol Evol. 22:1561–1568.
tive radiation in African cichlid fish. Nature. 513:375–381. Ho SY, Shapiro B, Phillips MJ, Cooper A, Drummond AJ. 2007. Evidence for 
Broad Institute. 2015. Picard version 1.141. [accessed 2015 December 16]. time dependency of molecular rate estimates. Syst Biol. 56:515–522.
Available from: https://broadinstitute.github.io/picard Hu X, Zhu W, Chen S, Liu Y, Sun Z, Geng T, Wang X, Gao B, Song C, Qin A, 
Callicrate  T, Dikow  R, Thomas  JW, Mullikin  JC, Jarvis  ED, Fleischer  RC; et al. 2016. Expression of the env gene from the avian endogenous retro-
NISC Comparative Sequencing Program. 2014. Genomic resources for the virus ALVE and regulation by miR-155. Arch Virol. 161:1623–1632.
endangered Hawaiian honeycreepers. BMC Genomics. 15:1098. IUCN. 2019. The IUCN Red List of Threatened Species version 2019-1. [ac-
Camacho  C, Coulouris  G, Avagyan  V, Ma  N, Papadopoulos  J, Bealer  K, cessed 2019 July 10]. Available from: https://www.iucnredlist.org
Madden  TL. 2008. BLAST+: architecture and applications. BMC Bio- James HF. 2004. The osteology and phylogeny of the Hawaiian finch radi-
informatics. 10:421. ation (Fringillidae: Drepanidini), including extinct taxa. Zool J Linn Soc. 
Carson HL, Clague DA. 1995. Geology and biogeography of the Hawaiian 141:207–255.
Islands. In: Wagner W, Funk V, editors. Hotspot archipelago. Washington James HF, Olson SL. 1991. Description of thirty-two new species of birds from 
(DC): Smithsonian Institution Press. p. 14–29. the Hawaiian Islands: part II. Passeriformes. Ornithological Monographs 
Cassin-Sackett L, Callicrate TE, Fleischer RC. 2019. Parallel evolution of gene 46. Washington (DC): American Ornithologists’ Union.
classes, but not genes: evidence from Hawai’ian honeycreeper populations Jarvis ED, Mirarab S, Aberer AJ, Li B, Houde P, Li C, Ho SY, Faircloth BC, 
exposed to avian malaria. Mol Ecol. 28:568–583. Nabholz B, Howard JT, et al. 2014. Whole-genome analyses resolve early 
Chaves JA, Cooper EA, Hendry AP, Podos J, De León LF, Raeymaekers JA, branches in the tree of life of modern birds. Science. 346:1320–1331.
MacMillan WO, Uy JA. 2016. Genomic variation at the tips of the adap- Jones  P, Binns D, Chang HY, Fraser M, Li W, McAnulla C, McWilliam H, 
tive radiation of Darwin’s finches. Mol Ecol. 25:5282–5295. Maslen J, Mitchell A, Nuka G, et al. 2014. InterProScan 5: genome-scale 
Cortes-Rodriguez N, Campana MG, Berry L, Faegre S, Derrickson SR, Ha RR, protein function classification. Bioinformatics. 30:1236–1240.
Dikow RB, Rutz C, Fleischer RC. 2019. Population genomics and struc- Joron M, Frezal L, Jones RT, Chamberlain NL, Lee SF, Haag CR, Whibley A, 
ture of the critically endangered Mariana Crow (Corvus kubaryi). Genes. Becuwe M, Baxter SW, Ferguson L, et al. 2011. Chromosomal rearrange-
10:187. ments maintain a polymorphic supergene controlling butterfly mimicry. 
Danecek  P, Auton  A, Abecasis  G, Albers  CA, Banks  E, DePristo  MA, Nature. 477:203–206.
Handsaker  RE, Lunter  G, Marth  GT, Sherry  ST, et  al.; 1000 Genomes Kahle KT, Kulkarni AV, Limbrick DD Jr, Warf BC. 2016. Hydrocephalus in 
Project Analysis Group. 2011. The variant call format and VCFtools. Bio- children. Lancet. 387:788–799.
informatics. 27:2156–2158. Kent WJ. 2002. BLAT–the BLAST-like alignment tool. Genome Res. 12:656–664.
del Hoyo J, Elliott A, Sargatal J, Christie DA, Kirwan G, editors. 2019. Hand- Knowlton  JL, Flaspohler  DJ, Rotzel  Mcinerney  NC, Fleischer  RC. 2014. First 
book of the birds of the world alive. Barcelona (Spain): Lynx Editions. record of hybridization in the Hawaiian honeycreepers: ‘i‘iwi (Vestiaria 
[accessed 2019 July 11]. Available from: http://www.hbw.com coccinea) × ‘apapane (Himatione sanguinea). Wilson J Ornithol. 126:562–
de Villemereuil P, Rutschmann A, Lee KD, Ewen JG, Brekke P, Santure AW. 568.
2019. Little adaptive potential in a threatened passerine bird. Curr Biol. Kozak K, McMillan WO, Joron M, Jiggins CD. 2018. Genome-wide admixture 
29:889–894.e3. is common across the Heliconius radiation. bioRxiv. doi:10.1101/414201
Dobrynin P, Liu S, Tamazian G, Xiong Z, Yurchenko AA, Krasheninnikova K, Lamichhaney S, Berglund J, Almén MS, Maqbool K, Grabherr M, Martinez-
Kliver S, Schmidt-Küntzel A, Koepfli KP, Johnson W, et al. 2016. Gen- Barrio A, Promerová M, Rubin CJ, Wang C, Zamani N, et al. 2015. Evo-
omic legacy of the African cheetah, Acinonyx jubatus. Genome Biol. lution of Darwin’s finches and their beaks revealed by genome sequencing. 
16:277. Nature. 518:371–375.
Eggert LS, Beadell JS, McClung A, McIntosh CE, Fleischer RC. 2009. Evolu- Lamichhaney S, Han F, Berglund J, Wang C, Almén MS, Webster MT, Grant BR, 
tion of microsatellite loci in the adaptive radiation of Hawaiian honey- Grant PR, Andersson L. 2016. A beak size locus in Darwin’s finches facili-
creepers. J Hered. 100:137–147. tated character displacement during a drought. Science. 352:470–474.
Elphick CS, Roberts DL, Reed JM. 2010. Estimated dates of recent extinctions Larkin  MA, Blackshields  G, Brown  NP, Chenna  R, McGettigan  PA, 
for North American and Hawaiian birds. Biol Conserv. 143:617–624. McWilliam  H, Valentin  F, Wallace  IM, Wilm  A, Lopez  R, et  al. 2007. 
Fleischer RC. 2009. Honeycreepers, Hawaiian. In: Gillespie R, Clague DA, edi- Clustal W and Clustal X version 2.0. Bioinformatics. 23:2947–2948.
tors. Encylopedia of islands. Berkeley (CA): University of California Press. Lechner  M, Findeiss  S, Steiner  L, Marz  M, Stadler  PF, Prohaska  SJ. 2011. 
p. 410-414 Proteinortho: detection of (co-)orthologs in large-scale analysis. BMC 
Fleischer RC, McIntosh CE, Tarr CL. 1998. Evolution on a volcanic conveyor Bioinformatics. 12:124.
belt: using phylogeographic reconstructions and K-Ar-based ages of the Lerner HR, Meyer M, James HF, Hofreiter M, Fleischer RC. 2011. Multilocus 
Hawaiian Islands to estimate molecular evolutionary rates. Mol Ecol. resolution of phylogeny and timescale in the extant adaptive radiation of 
7:533–545. Hawaiian honeycreepers. Curr Biol. 21:1838–1844.
Fleischer  RC, Tarr  CL, James  HF, Slikas  B, McIntosh  CE. 2001. Phylogen- Li H. 2013. Aligning sequence reads, clone sequences and assembly contigs 
etic placement of the po‘ouli Melamprosops phaeosoma based on mito- with BWA-MEM. arXiv. 1303.3997.
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022
32 Journal of Heredity, 2020, Vol. 111, No. 1
Li H, Durbin R. 2011. Inference of human population history from individual Smit AFA, Hubley R, Green P. 2013–2015. RepeatMasker open-4.0. [accessed 
whole-genome sequences. Nature. 475:493–496. 2017 October 13]. Available from: http://www.repeatmasker.org
Li  H, Handsaker  B, Wysoker  A, Fennell  T, Ruan  J, Homer  N, Marth  G, Song  L, Florea  L, Langmead  B. 2014. Lighter: fast and memory-efficient 
Abecasis G, Durbin R; 1000 Genome Project Data Processing Subgroup. sequencing error correction without counting. Genome Biol. 15:509.
2009. The sequence alignment/map format and SAMtools. Bioinformatics. Stamatakis A. 2014. RAxML version 8: a tool for phylogenetic analysis and 
25:2078–2079. post-analysis of large phylogenies. Bioinformatics. 30:1312–1313.
Li S, Li B, Cheng C, Xiong Z, Liu Q, Lai J, Carey HV, Zhang Q, Zheng H, Stanke  M, Diekhans  M, Baertsch  R, Haussler  D. 2008. Using native and 
Wei S, et al. 2014. Genomic signatures of near-extinction and rebirth of syntenically mapped cDNA alignments to improve de novo gene finding. 
the crested ibis and other endangered bird species. Genome Biol. 15:557. Bioinformatics. 24:637–644.
Lipman PW, Calvert AT. 2013. Modeling volcano growth on the Island of Ha- Stanke M, Schöffmann O, Morgenstern B, Waack S. 2006a. Gene prediction in 
waii: deep-water perspectives. Geosphere 9:1348–1383. eukaryotes with a generalized hidden Markov model that uses hints from 
Mack  KL, Nachman  MW. 2017. Gene regulation and speciation. Trends external sources. BMC Bioinformatics. 7:62.
Genet. 33:68–80. Stanke M, Tzvetkova A, Morgenstern B. 2006b. AUGUSTUS at EGASP: using 
Malinsky  M, Svardal  H, Tyers  AM, Miska  EA, Genner  MJ, Turner  GF, EST, protein and genomic alignments for improved gene prediction in the 
Durbin R. 2018. Whole-genome sequences of Malawi cichlids reveal mul- human genome. Genome Biol. 7(Suppl 1):S11.1–S11.8.
tiple radiations interconnected by gene flow. Nat Ecol Evol. 2:1940–1955. Supple MA, Hines HM, Dasmahapatra KK, Lewis JJ, Nielsen DM, Lavoie C, 
Marco-Sola  S, Sammeth  M, Guigó  R, Ribeca  P. 2012. The GEM mapper: Ray  DA, Salazar  C, McMillan  WO, Counterman  BA. 2013. Genomic 
fast, accurate and versatile alignment by filtration. Nat Methods. 9: architecture of adaptive color pattern divergence and convergence in 
1185–1188. Heliconius butterflies. Genome Res. 23:1248–1257.
Martin M. 2011. Cutadapt removes adapter sequences from high-throughput Sutton  JT, Helmkampf  M, Steiner  CC, Bellinger  MR, Korlach  J, Hall  R, 
sequencing reads. EMBnet J. 17:10-12. Baybayan P, Muehling J, GU J, Kingan S, et al. 2018. A high-quality, long-
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, read de novo genome assembly to aid conservation of Hawaii’s last re-
Garimella K, Altshuler D, Gabriel S, Daly M, et al. 2010. The genome ana- maining crow species. Genes. 9:393.
lysis toolkit: a mapreduce framework for analyzing next-generation DNA Tarr CL, Fleischer RC. 1993. Mitochondrial-DNA variation and evolutionary 
sequencing data. Genome Res. 20:1297–1303. relationships in the amakihi complex. Auk. 100:825–831.
Moyle RG, Oliveros CH, Andersen MJ, Hosner PA, Benz BW, Manthey JD, Toews  DP, Taylor  SA, Vallender  R, Brelsford  A, Butcher  BG, Messer  PW, 
Travers  SL, Brown RM, Faircloth BC. 2016. Tectonic collision and up- Lovette IJ. 2016. Plumage genes and little else distinguish the genomes of 
lift of Wallacea triggered the global songbird radiation. Nat Commun. hybridizing warblers. Curr Biol. 26:2313–2318.
7:12709. Tollis M, Hutchins ED, Stapley J, Rupp SM, Eckalbar WL, Maayan I, Lasku E, 
Pratt HD. 2005. The Hawaiian honeycreepers: Drepanidinae. Oxford (UK): Infante CR, Dennis SR, Robertson JA, et al. 2018. Comparative genomics 
Oxford University Press. reveals accelerated evolution in conserved pathways during the diversifi-
Ralph CJ, Fancy SG. 1996. Aspects of the life history and foraging ecology of cation of anole lizards. Genome Biol Evol. 10:489–506.
the endangered akiapolaau. Condor. 98:312–321. Treangen  TJ, Salzberg  SL. 2012. Repetitive DNA and next-generation 
Reed  RD, Papa  R, Martin  A, Hines  HM, Counterman  BA, Pardo-Diaz  C, sequencing: computational challenges and solutions. Nat Rev Genet. 
Jiggins CD, Chamberlain NL, Kronforst MR, Chen R, et al. 2011. optix 13:36–46.
drives the repeated convergent evolution of butterfly wing pattern mim- UniProt Consortium. 2019. UniProt: a worldwide hub of protein knowledge. 
icry. Science. 333:1137–1141. Nucleic Acids Res. 47:D506–D515.
Roberts DL, Elphick CS, Reed  JM. 2009. Identifying anomalous reports of van Riper C III, van Riper SG, Goff ML, Laird, M. 1986. The epizootiology 
putatively extinct species and why it matters. Conserv Biol. 24:189–196. and ecological significance of malaria in Hawaiian (USA) land birds. Ecol 
Salzburger W. 2018. Understanding explosive diversification through cichlid Monogr. 56:327–344.
fish genomics. Nat Rev Genet. 19:705–717. Warner RE. 1968. The role of introduced diseases in the extinction of the en-
Samuel  MD, Hobbelen  PHF, DeCastro  F, Ahumada  JA, LaPointe  DA, At- demic Hawaiian avifauna. Condor. 7: 101–120.
kinson  CT, Woodworth  BL, Hart  PJ, Duffy  DC. 2011. The dynamics, Warren  WC, Clayton  DF, Ellegren  H, Arnold  AP, Hillier  LW, Künstner  A, 
transmission, and population impacts of avian malaria in native Hawaiian Searle S, White S, Vilella AJ, Fairley S, et al. 2010. The genome of a song-
birds: a modeling approach. Ecol Appl. 21:2960–2973. bird. Nature. 464:757–762.
Samuel MD, Woodworth BL, Atkinson CT, Hart PJ, LaPointe DA. 2015. Avian Waterhouse RM, Seppey M, Simão FA, Manni M, Ioannidis P, Klioutchnikov G, 
malaria in Hawaiian forest birds: infection and population impacts across Kriventseva EV, Zdobnov EM. 2017. BUSCO applications from quality 
species and elevations. Ecosphere. 6:1–21. assessments to gene prediction and phylogenomics. Mol Biol Evol. 
Simpson JT, Wong K, Jackman SD, Schein JE, Jones SJ, Birol I. 2009. ABySS: a 35:543–548.
parallel assembler for short read sequence data. Genome Res. 19:1117–1123. Wittkopp  PJ, Kalay  G. 2012. Cis-regulatory elements: molecular mechan-
Slater GS, Birney E. 2005. Automated generation of heuristics for biological isms and evolutionary processes underlying divergence. Nat Rev Genet. 
sequence comparison. BMC Bioinformatics. 6:31. 13:59–69.
Smeds  L, Qvarnström  A, Ellegren  H. 2016. Direct estimate of the rate of Wolf JB, Ellegren H. 2017. Making sense of genomic islands of differentiation 
germline mutation in a bird. Genome Res. 26:1211–1218. in light of speciation. Nat Rev Genet. 18:87–100.
Smeds L, Warmuth V, Bolivar P, Uebbing S, Burri R, Suh A, Nater A, Bureš S, Zuccon D, Prŷs-Jones R, Rasmussen PC, Ericson PG. 2012. The phylogenetic 
Garamszegi  LZ, Hogner  S, et  al. 2015. Evolutionary analysis of the relationships and generic limits of finches (Fringillidae). Mol Phylogenet 
female-specific avian W chromosome. Nat Commun. 6:7330. Evol. 62:581–596.
Downloaded from https://academic.oup.com/jhered/article/111/1/21/5625504 by smithsonia3 user on 31 May 2022