Mitochondrial Genomes from New Zealand’s Extinct Adzebills (Aves: Aptornithidae: Aptornis ) Support a Sister-Taxon Relationship with the Afro-Madagascan Sarothruridae

: The recently extinct New Zealand adzebills (Aptornithidae, Aptornis spp.) were an enigmatic group of large flightless birds that have long eluded precise taxonomic assignment as they do not closely resemble any extant birds. Adzebills were nearly wingless, weighed approximately 16–19 kg, and possessed massive adze-like reinforced bills whose function remains unknown. Using hybridisation enrichment and high-throughput sequencing of DNA extracted from subfossil bone and eggshell, near-complete mitochondrial genomes were successfully assembled from the two Quaternary adzebill species: the North Island Adzebill ( Aptornis otidiformis ) and South Island Adzebill ( A. defossor ). Molecular phylogenetic analyses conﬁrm that adzebills are members of the Ralloidea (rails and allies) and are sister-taxon to the Sarothruridae, which our results suggest comprises the Madagascan wood rails ( Mentocrex , two likely sp.) in addition to the tiny (<50 gram) rail-like Afro-Madagascan ﬂufftails ( Sarothrura , 9 spp.). Node age estimates indicate that the split between adzebills and Sarothruridae occurred ~39.6 Ma, suggesting that the ancestors of the adzebills arrived in New Zealand by long-distance dispersal rather than continental vicariance. This newly identiﬁed biogeographic link between physically distant New Zealand and Afro-Madagascar, echoed by the relationship between the New Zealand kiwi (Apterygiformes) and Madagascan elephant-birds (Aepyornithiformes), suggests that the adzebill’s near relatives were formerly more widespread. Quaternary adzebill species land-bridge and are the result of a relatively recent dispersal from the South Island, from which the earliest (Miocene) adzebill fossil has been described.


Introduction
The islands of New Zealand (NZ) are exceptional for their long oceanic isolation from the supercontinent Gondwana (>52 Ma) [1], large size (>260,000 km 2 ), and absence of mammals (except bats) among the recent native land-fauna. NZ's native biodiversity is rich in phylogenetically distinct lineages, and extraordinary examples of gigantism and flightlessness among the avifauna. Important examples include the nine species of megaherbivore ratite moa, weighing up to 250 kg (Dinornithiformes) [2,3], four of the world's five known flightless songbirds [4], the world's heaviest parrot (also flightless) [5], and the largest-known eagle (Aquila moorei) [6]. However, widespread extinction (>40% of endemic bird species) [3,7] following human settlement in NZ around 1270 AD [8] means that phylogenetic studies on many NZ taxa have been dependent on advances in ancient DNA (aDNA) research [9]. As a result, the evolutionary histories of many of NZ's unique animal taxa, including the adzebill (Aptornis: Aptornithidae), remain unresolved or inconclusive.
Adzebills were large, specialised, flightless birds endemic to NZ, which have been a taxonomic and ecological enigma since they were first described by Richard Owen in 1844 [10][11][12][13]. The genus Aptornis includes two recently extinct species: the North Island Adzebill (A. otidiformis), typically reaching around 16 kg; and the larger South Island Adzebill (A. defossor), typically reaching around 19 kg (though a maximum size of 25 kg has been suggested) [2]. In addition, fragments of a fossil adzebill (A. proasciarostratus), strikingly similar to recent adzebills, have also been described from the Early Miocene (16)(17)(18)(19) Ma) deposits near St Bathans in the southern South Island, suggesting a long history of the group in NZ [14]. Both recent species were nearly wingless and had a robust and specialised morphology. The massive beak and skull, as well as the neck vertebrae, were heavily reinforced and likely supported powerful muscles. Adzebills also had robust legs and feet, which may have been used with the massive bill to excavate or immobilize food or prey. Stable isotopes from adzebill bone specimens confirm that adzebills had a high trophic niche and were likely predators or scavengers, although their exact feeding strategy remains unknown [2,15]. Archaeological deposits confirm that adzebills were hunted by early Maori [16][17][18], who also rapidly cleared the dry, lowland podocarp forests in eastern NZ [19] which were the birds' main habitat during the Holocene [2,15]. Both activities led to the extinction of the recent adzebill species around the same time as the moa, circa 1500 CE [2,20].
The highly derived morphology of adzebills has long complicated their classification, though they have usually been placed in Gruiformes along with several other bird families found in former fragments of Gondwana (especially New Caledonia, Madagascar, and South America). It has often been proposed that the ancestor of these "gruiform" bird families, including adzebills, inhabited Gondwana prior to its break up and that their present biogeography reflects continental vicariance [21]. Specifically, morphological analyses have suggested adzebills are most closely related to the New Caledonian Kagu (Rhynochetos jubatus) [21][22][23][24], consistent with the possible existence of emergent land connecting New Zealand and New Caledonia during the Paleogene [25,26]. However, recent phylogenetic studies have found that most traditional gruiforms do not form a monophyletic group, except for two superfamilies (each comprising three extant families): the Gruoidea (cranes and allies) and Ralloidea (rails and allies) [27][28][29]. For example, the Kagu and the South American Sunbittern Eurypyga helias (together comprising Eurypygiformes) appear to form the sister-taxon to tropicbirds (Phaethontiformes) [24,29,30]. Previous analyses of a 673 bp fragment of the mitochondrial 12S rRNA gene from the South Island Adzebill suggested that adzebills were members of the Ralloidea in the Gruiformes and unrelated to kagu [31], although insufficient data and limited taxon sampling prevented confident identification of the adzebills' precise phylogenetic affinities. These alternative and incompatible hypotheses for the classification of adzebills-as members of either Eurypygiformes or Gruiformes-have widely disparate implications for their biogeographic and temporal origins.
If adzebills are ralloids, and not eurypygiforms, then the question becomes: "to what living ralloids are they most closely related?" Any attempts to answer this question require a robust phylogenetic framework for Ralloidea, which historically contains an uncertain number of families (with varying taxonomic contents). For example, the flufftails (Sarothrura spp.) were previously considered members of the ralloid family Rallidae (the rails) but are now recognised as comprising their own family (Sarothruridae) [32] more closely related to the finfoots (Heliornithidae) [24,28]. Similarly, the Madagascan wood rails (Mentocrex kioloides sspp.), which are presently classified as Rallidae [32], were recently suggested to be closer relatives of Sarothruridae [33].
In this study, we present the first near-complete mitochondrial genome sequences from both recent species of adzebill. We also present new genetic data from several key gruiform lineages of uncertain affinity, including members of the Rallidae. We analyse these new sequences alongside all available mitochondrial data for ralloids, as well as the kagu, to confidently resolve the phylogenetic position of the adzebills and estimate the timeline for their evolution.

Samples
We sampled femora from one North Island (NI) Adzebill (A. otidiformis) and one South Island (SI) Adzebill (A. defossor), both held in the collections of the Museum of New Zealand Te Papa Tongarewa (NMNZ). The North Island specimen (NMNZ DM8046) was collected from a cavern near Coonoor, Wairarapa, south-eastern North Island (40.4 • S, 176.1 • E). Worthy [34] inferred that the deposits in this cavern had a maximum age of "a few thousands of years at most". The SI specimen (NMNZ S23033) was collected from the Honeycomb Hill Cave System, north-western South Island (41.1 • S, 172.2 • E). Adzebills were most common in Pleistocene deposits in the north-western South Island (including Honeycomb Hill system) [2,15,35,36], so S23033 may be of Pleistocene age (i.e., >12 ka).
In addition to the two bone samples, we also analysed three putative adzebill eggshell fragments collected from a rockshelter near Hukanui, central North Island (39.3 • S, 176.5 • E), held in the collections of the Canterbury Museum. The collection (CMC AV18536) consists of about 40 shards of eggshell, 20 of which were identified as having the same curvature and pore structure as the fragments confirmed as Aptornis (see below). These Aptornis-like fragments varied between 3 and 8 mM in diameter and averaged between 610 and 650 µM in thickness. The fragments were off-white in colour and displayed a low density of round pores on the surface. Bone and sediment samples from this rockshelter [37] have been shown using AMS (Accelerator Mass Spectrometry) and OSL (Optically-Stimulated Luminescence) dating to post-date the Taupo Ignimbrite (c. 1850 yr BP).
The DNA extraction and library preparation for the subfossil bone and eggshell samples were performed in purpose-built, physically isolated, ancient DNA laboratories. The bone samples were processed at the Australian Centre for Ancient DNA (ACAD), University of Adelaide, while the eggshell fragments were processed at the Ancient DNA Laboratory, Murdoch University.
We also sampled a number of extant gruiforms to serve as outgroups in our phylogenetic analyses (Table S1). Muscle of the extant Gray-winged Trumpeter, Psophia crepitans (Gruoidea: Psophiidae), was sampled from a deceased individual from the National Zoological Park, USA (species assignment was confirmed by comparison to the available P. crepitans CytB sequence, GenBank accession D8Q4900, which yielded a 99.96% pairwise match). In addition, we also sourced skin samples from museum specimens of the Limpkin Aramus guarauna (Gruoidea: Aramidae), the Red-chested Flufftail Sarothrura rufa (Ralloidea: Sarothruridae), and the Grey-throated Rail Canirallus oculeus (Ralloidea: Rallidae) (Table S1). Although no genetic data from C. oculeus has previously been analysed, it is often considered congeneric with Mentocrex, and may therefore be a close relative of Sarothruridae. Samples from extant species were processed in a clean-room at the University of Adelaide or the South Australian Regional Facility for Molecular Ecology and Evolution (SARFMEE) molecular biology labs.

Bone Samples
To control for contaminants, the exterior surfaces (~1 mm) of both adzebill bone samples were removed using a Dremel tool. Bone samples were powdered using a mikro-dismembrator (Sartorius). Around 0.2-0.3 g of bone powder from each sample was digested by rotational incubation at 37 • C overnight in 4 mL of 0.5M EDTA pH 8.0 followed by a second round of incubation at 55 • C with 60 µL of proteinase-K. The resulting solution was bound to silica, washed, and eluted using the silica-based method of Brotherton et al. [38].

Tissue Samples
DNA was isolated from samples of extant species using a DNeasy Blood & Tissue kit (for skin) (QIAGEN, Hilden, Germany) or a modified salting-down extraction (for muscle), where~0.3 g tissue was mixed with 300 µL cell lysis solution, 3 µL proteinase K and incubated at 55 • C [39]. The high molecular weight DNA extracted from the grey-winged trumpeter sample was sheared to~300 bp using a focused-ultrasonicator (Covaris). DNA extracted from the other extant species was not sheared, as the specimens were museum skins and the DNA was consequently already fragmented.

Eggshell Samples
Exposed surfaces of the eggshell fragments were abraded with a Dremel tool at low speed and the resulting powder was discarded to minimise contamination. The cleaned eggshell fragments (each 50-100 mm 2 in size) were then further abraded with the Dremel to produce additional powder, which was extracted using the digestion and silica-binding method of Oskam et al. [40].

Bone and Tissue Samples
Extracts were blunt-end repaired, and custom adapters ligated, using the library preparation protocol of Meyer and Kircher [34]. Each adapter (5' and 3') contained a unique 7 mer index, allowing both for identification of amplified DNA and downstream removal of contaminant sequencing reads. Following library preparation, each sample was amplified using PCR in eight separate reactions to reduce amplification bias. Each reaction contained 1 × PCR buffer, 2.5 mM MgCl2, 1 mM dNTPs, 0.5 mM primer, 1.25 U AmpliTaq Gold, 2 µL DNA extract. Reactions were subjected to the following thermocycling regime: 94 • C 12 min; 13 of 94 • C for 30 s, 60 • C for 30 s, 72 • C for 40 s (plus an additional 2 s per cycle); and a final extension of 72 • C for 10 min. Individual PCR products for each library were pooled following amplification, and purified using AMPure magnetic beads (Agencourt).
All libraries were enriched for host mitochondrial DNA using commercially synthesised biotinylated 80 mer RNA baits (Arbor Biosciences, Ann Arbor, MI, USA). Baits were designed using a wide range of published whole mitochondrial genome sequences (not including control region/D-loop) of a variety of avian taxa, including palaeognaths, galloanseres, and neoavians [41][42][43][44]. DNA-RNA hybridisation enrichment was performed using 200 ng of library following the manufacturer's recommended protocol (myBaits v1), with the exception that the incubation was extended to 44 h (3 h at 60 • C, 12 h at 55 • C, 12 h at 50 • C, 17 h at 55 • C). Following incubation, baits were immobilised on magnetic MyOne Streptavidin Beads (Life Technologies, Christchurch, New Zealand). The baits were washed once with 1 × SCC and 0.1% SDS (15 min at room temperature), and twice with 0.1 × SCC and 0.1% SDS (10 min at 50 • C), then resuspended in 0.1 M NaOH pH 13.0. The resulting enriched library was purified using a Minelute spin-column (QIAGEN), and subjected to a further round of PCR (12 cycles, eight reactions, using above recipe). Enriched libraries were subjected to a final round of PCR (seven cycles, five reactions, using above recipe), using fusion primers to add full-length Illumina sequencing adapters for sequencing [45]. Libraries were diluted to 2 nM (quantified using KAPA qPCR) and sequenced on an Illumina MiSeq using 2 × 150 (300 cycle, paired-end) sequencing chemistry.

Eggshell Samples
A fragment of the mitochondrial 12S rRNA gene was amplified and sequenced from the eggshell DNA extractions using the primers 12Sa and 12Sh [46], which confirmed the identity of two samples as A. otidiformis. The third eggshell DNA extract was identified as Little Bush Moa (Anomalopteryx didiformis) and was excluded from further analyses. The eggshell libraries were prepared for sequencing using a Roche 454 GS Junior platform. Roche RL017 adapters were ligated to pooled DNA from the two eggshell extracts that were positively identified as A. otidiformis using the New England Biosciences Next Quick DNA sample prep reagent set 2 (Protocol 1). Unligated adapter oligos were removed by undertaking the small fragment removal (E6116) twice.
RNA baits (Arbor Biosciences) were produced in 100 mer lengths with 50 bp tiling density based on published mitochondrial genomes from 10 avian taxa (Table S2), excluding the control region/D-loop. A second set of RNA baits were produced based on nuclear genes from six gruiform species published by Hackett et al [28]. An equimolar mixture of these two baits sets was made and used to enrich our library, following the manufacturer's protocol (myBaits v1) with the exception that the library was eluted from the MyOne beads using a Tris-based buffer rather than NaOH. Small fragments were removed from the resulting captured library by performing two sequential rounds of clean-up with Agencourt AMPure XP beads, before eluting a final volume of 40 µL.
Ion Torrent (ThermoFisher, Waltham, MA, USA) A and B adapters were ligated to the enriched library during a 25-cycle amplification by incorporating Ion Torrent adapters into the 454 primers. As a result, the post amplification library contained sequences for both Ion Torrent and 454 adapters with an Ion Torrent keypass between the two. Sequencing was performed on an Ion Torrent Personal Genome Machine (PGM) with 316 chip array in accordance with the manufacturer's protocol.

Bone and Tissue Samples
Reads belonging to each library were demultiplexed according to their unique combination of 5' and 3' barcode sequences using Sabre v1.0 (https://github.com/najoshi/sabre) (default parameters; no mismatches allowed). Adapter sequences were removed and paired-end reads were collapsed (where the forward and reverse reads overlapped) using AdapterRemoval v2.1.2 [47]. Low quality bases were trimmed (<Phred20 -minquality 4) and collapsed reads shorter than 25 bp were discarded (-minlength 25). Read quality was visualised using fastQC v0.10.1 (http://www.bioinformatics.bbsrc. ac.uk/projects/fastqc) before and after trimming to make sure the trimming of adapters was efficient. Only the collapsed reads from the bone and tissue samples were used in downstream analyses.

Eggshell Samples
Ion Torrent PGM post-run analysis identified and removed the Ion Torrent keypass and adapter sequences. We then used Geneious v8.1.9 to remove the 454 adapter sequences by breaking the 3.2 million reads into 500,000 read batches (due to memory limitations) and searching for a 25 bp motif corresponding to the 454 adapter sequence. In addition, we trimmed the last 35 bp of all reads to ensure that any residual adapter was removed from the opposite end of the molecule.

Genome Assembly
To create consensus sequences for outgroup species and our highest-quality adzebill specimen (DM8046), we iteratively mapped collapsed reads from each species to the published NZ Weka (Gallirallus australis) mitochondrial genome (KF701060) using the Geneious Read Mapper in Geneious v9.1.6 [48,49] with the following parameter values: maximum mismatches per read = 15%, minimum mapping quality = Phred 30, maximum gaps per read = 10% read length, maximum gap length = 10 bp, minimum overlap = 25 bp, word size = 8. We continued mapping iteratively using these parameter values until no additional reads were obtained and then created a majority-rule consensus sequence after checking the alignment by eye. We then used BWA [50] (aln -l 1024, -o 2) to perform a more stringent round of mapping for each species using the majority-rule consensuses as reference sequences. Mapping stringency in BWA was relaxed slightly for adzebills versus the extant outgroups (−n 0.01 versus −n 0.04) to account for the presence of miscoding lesions in our aDNA data. Final 75% majority consensus sequences were then generated in Geneious after removing duplicate reads using 'FilterUniqueSAMCons.py' [51]. We only called bases for sites covered by ≥3 reads (appropriate IUPAC degenerate bases were called instead where coverage was insufficient or no >75% majority base was observed). To generate consensus sequences for our other two adzebill specimens (S23033 and AV18536), we mapped the collapsed reads from S23033 and all reads from the eggshell library to the majority-rule consensus for DM8046 using BWA (aln -l 1024, -o 2, -n 0.01). Final consensus sequences were generated for these two specimens as above, with the exception that "MarkDuplicates.py" was used to deduplicate the reads from the eggshell library, as the Ion Torrent PGM data was single-end only (so reads derived from the same template molecule may be of different lengths). Details of all final consensus sequences are provided in Table S1. We used MapDamage v2.0.6 [52] to ensure that damage patterns in the adzebill data were consistent with authentic ancient DNA (Figures S1-S3).

Taxon Sampling and Sequence Alignment
We made three sequence alignments with varying taxon-sampling and levels of data "completeness" (see Table S3 and description below for sample accession numbers and genes selected). All alignments were of identical length (15,140 bp) and differed only in respect to taxa included. In these alignments, we used only the most complete sequence from our two North Island Adzebill samples (DM8046), as the two sequences were nearly identical for alignment columns where both were unambiguously called. For other taxa where the same locus had been sequenced for multiple individuals we created and used a 90% majority rule consensus. Exceptions were made for several taxa, such as the New Guinea Flightless Rail (Megacrex inepta), where pairwise comparison revealed high levels of divergence between sequences supposedly from individuals belonging to the same species (in these cases sequences from different specimens/studies were considered independently). Alignments were performed using MAFFT v7.017 [53,54] (default settings) as implemented in Geneious, and sequences were partitioned into protein-coding genes, tRNAs, and rRNAs. We checked the alignments by eye and removed poorly aligning regions around large insertions/deletions. Only one taxon in our alignments, the Black-browed Albatross Thalassarche melanophris (Procellariiformes), had an observable gene duplication (two adjacent copies of the region comprising ND6, control region and nearby tRNAs): a pattern widespread in the genus Thalassarche [55]. In the case of T. melanophris, only the duplicated gene region next to Cytochrome B was retained in our alignments.
Our first alignment (alignment A), comprised 12 complete mitochondrial genome sequences: the North Island adzebill; the Kagu (Rhynochetos jubatus); one member of each gruiform family, including the Eurasian Coot Our second alignment (alignment B) was focused on Gruiformes and comprised mostly complete mitochondrial genome sequences (excluding the control region), together representing all the major gruiform clades identified by previous phylogenetic studies [27,33,56] (Table S3). We also included the Southern Rockhopper Penguin and Black-browed Albatross as outgroups. As Krajewski et al. [57] have already comprehensively resolved the mitochondrial phylogeny of Gruidae, in our study we elected to include only the Black-crowned Crane (Balearica pavonina) to represent the Balearicinae and the Common Crane and Siberian Crane (Leucogeranus leucogeranus, the basal gruine species) to represent the Gruinae. Although not represented by complete mitochondrial genomes, we included all available data from the African Finfoot and Masked Finfoot (Podica senegalensis and Heliopais personatus, respectively) of the Heliornithidae and the Madagascan wood rails (Mentocrex sspp.) of the Rallidae, so that all ralloid families were represented.
Our final alignment (alignment C) comprised the exact same taxa and data from alignment B but with the addition of more taxa for which only one or several mitochondrial genes sequences were available. Thus, of our three alignments, alignment C had the most complete taxon sampling but also the most missing data. Overall, alignment C covered all gruiform families, including both Gruidae subfamilies (Gruinae, 2/3 genera, 2/3 species; Balearicinae, 1/1 genera, 1/2 species), Psophiidae (1/1 genera, 1/3 species), Aramidae (1/1 genera, 1/1 species), Heliornithidae (3/3 genera, 3/3 species) and Sarothruridae (1/1 genera, 2/9 species). The taxonomy of the speciose family Rallidae is in continual flux and so, for reference, in the present study we follow the classification of Dickinson and Remsen [32], which recognises 146 species in 40 genera (extant and recently extinct). Of the taxa recognised by Dickinson and Remsen [32], alignment C includes representatives from 98 species and 32 genera.

Molecular Phylogenetic Analyses
To test for saturation of 3rd codon positions of protein coding genes, we used DAMBE v7.0.28 [58,59]. Using alignment B, but omitting Mentocrex, Heliopais, and Podica (which lacked complete mitochondrial genomes), we partitioned our dataset into 1st, 2nd, and 3rd codon positions of protein-coding genes, and into RNA-coding sequences. The saturation test by Xia et al. [60] found that 3rd codon positions were significantly more saturated than all other partitions (and were the only positions for which significant saturation was detected), although they still contained phylogenetic information (p < 0.01, Table S4). Saturation plots of transitions/transversions versus divergence of all partitions were generated using an F84 substitution model ( Figure S4). 1st codon, 2nd codon, and RNA-coding sequences revealed linear plots for both transitions and transversions typical of unsaturated sites [60]. However, the saturation plot of 3rd codon positions demonstrated a linear model for transversions only, whereas transitions quickly became saturated with increasing divergence. As a result, all 3rd codon position nucleotides in our alignments were converted into the IUPAC ambiguity codes R (purines) and Y (pyrimidines), thus removing transitions at 3rd codon positions from all subsequent analyses.
To identify appropriate partition schemes and substitution models for each of our three alignments we used PartitionFinder v2.1.1 [61,62]. We provided PartitionFinder with 51 putative partitions, comprising 1st and 2nd codon positions of individual protein-coding genes, plus the individual rRNA tRNAs coding genes. 3rd codon positions (RY-coded) were combined into a single partition. We used the greedy algorithm and chose the best-supported scheme according to the Bayesian Information Criterion (Table S5).
Phylogenetic trees were estimated using Bayesian estimates as implemented through MrBayes v3.2.6 [63,64], using the partition schemes and substitution models suggested by PartitionFinder. These used a total of 1 × 10 6 (alignment A and B) and 5 × 10 6 generations (alignment C). A total of eight separate chains were used, and convergence between all chains was assessed and confirmed using the standard convergence diagnostic (standard deviation of split frequencies). The first 25% of sampled trees were discarded as a burn-in to estimate tree topology.
In addition, we performed a maximum likelihood phylogenetic analysis using RAxML v8.2.0 [65] on alignment A only (the reduced taxon set which included the Kagu Rhynochetos jubatus). Our RAxML analysis comprised a maximum likelihood search for the best scoring tree from 1000 bootstrap replicates (-f a). Subsequently, we used RAxML to perform a Shimodaira-Hasegawa test [66,67] to compare our best tree topology to several alternative topologies (-f h): (1) a tree where Aptornis otidiformis was sister taxon to R. jubatus; (2) a tree where A. otidiformis was sister-taxon to Psophia crepitans; (3) a tree where A. otidiformis was sister-taxon to a clade comprising P. crepitans, Grus grus, and Aramus guarauna; and (4) a tree where Aptornis was sister-taxon to Fulica atra.

Morphological Phylogenetic Analyses
In addition to our molecular analyses, we also reanalysed previously published morphological character matrices [22,23]. Livezey and Zusi's [23] results suggested a sister-relationship between Aptornis and Rhynochetos, which conflicted with the results of our molecular analyses (which supported Aptornis as a member of Gruiformes; see Results section). Consequently, we sought to identify any morphological characters that supported our molecular results. We reduced Livezey and Zusi's dataset of 2954 characters from 188 taxa down to 456 parsimony informative characters for 10 taxa: Diomedea (Procellariiformes), Spheniscus (Sphenisciformes), Eurypyga (Eurypygidae), Rhynochetos, Psophia, Aramus, Grus, Porphyrio (Rallidae), Heliornis, and Aptornis. Using accelerated transformation in PAUP 4.0a [68] we mapped the characters onto constraint trees (following [29,30]) in which Aptornis was placed as sister taxon to its nearest sampled relative, Heliornis (as identified by our molecular analysis).
We attempted to identify any additional unrecognised morphological apomorphies diagnosing the phylogenetic position of the adzebills (as revealed by our molecular analyses) by reanalysing Livezey's [22] morphological character matrix. Livezey [22] compiled an extensive dataset comprising 225 taxa sampled for 381 primarily osteological characters and 189 integument characters. Our reanalysis was performed using PAUP 4.0a [68], and comprised parsimony searches employing 100 random additions, the default heuristic search settings, and 1000 bootstrap replicates. The morphological dataset was analysed using the "backbone constraints" command, constraining the relationships among the 15 taxa represented in our molecular dataset to be consistent with the topology of the trees in Figures 1-3. Other taxa were left unconstrained and their phylogenetic positions were inferred from the morphological data. The "describe trees/apolist" command in PAUP was used to find ambiguous and unambiguous apomorphies defining the clade containing Aptornithidae and its nearest relatives in the bootstrap tree.

Fossil Calibrations and Molecular Clock Analyses
Molecular dating analyses were performed on alignment B using the BEAST software package v1.8.4 (which includes the complementary programs BEAST, LogCombiner, and TreeAnnotater) [69]. The analysis was run four times independently with BEAST v1.8.4 using a birth-death tree prior and a single relaxed lognormal clock model. Each chain was run for 2 × 10 8 generations, sampling every 2 × 10 4 generations, with the first 10% of samples discarded as burn-in. The four BEAST runs used different starting trees, each of which had the same topology as the consensus tree obtained from the MrBayes analysis (but had varying branch lengths). We used the partitions and substitution models suggested by PartitionFinder (Table S4). The trees from each run were combined using LogCombiner v1.8.4 and summarised in TreeAnnotater v1.8.4 to create a maximum clade credibility tree.  We calibrated our molecular dating analyses by constraining the age of four nodes in accordance with the fossil record: 1. We followed Claramunt and Cracraft [24] in constraining the split between the Sphenisciformes and Procellariiformes according to an exponential distribution with a mean of 11 Ma and offset of 61 Ma. The offset (minimum bound) was determined by the age of the oldest unequivocal member of the clade, the fossil penguin Waimanu manneringi [70] from the New Zealand Waipara Greensand. Claramunt and Cracraft [24] determined the mean value of the exponential distribution (11 Ma) empirically using the ages of the oldest fossil representatives of the clade (n = 6) from each continent (assuming these fossils were drawn from a uniform temporal distribution).
2. We also followed Claramunt and Cracraft [24] in constraining the split between Gruoidea and Ralloidea according to an exponential distribution with a mean of 8.5 Ma and offset of 52 Ma. Here the offset was determined by Messelornis nearctica [71], the oldest unequivocal member of the Ralloidea. As above, Claramunt and Cracraft [24] determined the value of the mean (8.5 Ma) based on the temporal distribution of the oldest members of the clade (n = 7) from each continent.
3. The well-preserved fossil crane Balearica exigua (Balearicinae) dated to >11.93 Ma [72], as well as similarly-aged subfossil cranes assigned to the genus Palaeogrus (Gruinae) [73], provide a conservative minimum age for the crane subfamilies Balearicinae and Gruinae. Consequently, we constrained the divergence between Balearicinae and Gruinae according to a uniform distribution with a minimum of 12 Ma and a maximum of 65.5 Ma. The maximum bound for this node assumes that all neoavian families and subfamilies diverged during the Tertiary, supported by the lack of crown members of neoavian families prior to the Paleocene and consistent with the results of most recent molecular studies.
4. The fossil adzebill Aptornis proasciarostratus from the New Zealand St Bathans Fauna (16-19 Ma) provides a robust minimum bound for the divergence between the Aptornis lineage and its nearest living relatives [14]. As for constraint 3 above, the beginning of the Tertiary provides a conservative upper bound for family level divergences within Neoaves. Consequently, we constrained the stem-age of Aptornis according to a uniform distribution with a minimum of 16 Ma and a maximum of 65.5 Ma.   Phylogenetic analyses of alignment C, which included more taxa but more missing data than alignments A and B, confirmed with high support that all sequenced Rallidae (sensu Dickinson and Remsen [32], except Mentocrex) comprised a monophyletic group that was sister taxon to the heliornithoid clade. However, missing data meant that the phylogenetic position of many species was poorly resolved (Figure 3). For example, the Chestnut-headed Crake Rufirallus castaneiceps was represented only by a short and relatively conserved fragment of 12S rRNA not covered by its most likely close relatives (i.e. R. viridis or Micropygia schomburgkii), which prevented us from confidently resolving the phylogenetic position of this species within Rallidae (as a result, we omitted R.

Mitochondrial Tree Topology
Our Bayesian phylogenetic analyses of alignments A and B, which included predominantly complete mitochondrial genome sequences, resulted in extremely well resolved trees (Figures 1  and 2; all nodes with Bayesian posterior probability of 1.0 unless stated otherwise). We confirmed that adzebills (Aptornithidae) were crown gruiforms, and more-specifically found that Aptornis was sister-taxon to a clade comprising Mentocrex and Sarothrura. Alternative phylogenetic positions for the adzebill, including a tree in which Aptornis and the kagu (Rhynochetos jubatus) were sister-taxa, could be rejected at p < 0.01 by a Shimodaira-Hasegawa test applied to alignment A (Figure 1). Our analyses further suggested that Aptornis, Mentocrex, and Sarothrura formed a clade that was sister to Heliornithidae. For brevity, we hereafter refer to members of the Aptornis-Mentocrex-Sarothrura-Heliornithidae clade (the sister clade to Rallidae) as "heliornithoids" since Heliornis has nomenclatural priority.
Placement of Mentocrex among the heliornithoids, as previously suggested by García-R et al. [33], is incompatible with its Rallidae classification, and we therefore recommend Mentocrex be transferred to the Sarothruridae. A previous study of morphological data suggested a close relationship (even congeneric status) between Mentocrex and Canirallus [22], suggesting that Canirallus may also not be in the Rallidae. However, our new data from the African Grey-throated Rail Canirallus oculeus allow us to refute this relationship. We reveal that Canirallus is an early-diverging member of Rallidae (Figures 2 and 3) and thus should remain classified as a rallid. Otherwise, the results of our phylogenetic analyses based on alignments A and B generally recapitulate previous studies [27,28,33,55,74,75]. In particular, our new mitochondrial genomes from Psophia and Aramus provided improved resolution for relationships among gruoids, which we confirmed to be monophyletic and sister-taxon to Ralloidea. Our phylogenetic analyses further confirmed that the basal gruoid family is Psophiidae (Bayesian posterior probability = 0.95), with Gruidae and Aramidae as sister-taxa, in agreement with the results of other molecular studies [27,28,74].
Phylogenetic analyses of alignment C, which included more taxa but more missing data than alignments A and B, confirmed with high support that all sequenced Rallidae (sensu Dickinson and Remsen [32], except Mentocrex) comprised a monophyletic group that was sister taxon to the heliornithoid clade. However, missing data meant that the phylogenetic position of many species was poorly resolved (Figure 3). For example, the Chestnut-headed Crake Rufirallus castaneiceps was represented only by a short and relatively conserved fragment of 12S rRNA not covered by its most likely close relatives (i.e., R. viridis or Micropygia schomburgkii), which prevented us from confidently resolving the phylogenetic position of this species within Rallidae (as a result, we omitted R. castaneiceps from our final analysis). In addition, our analyses of alignment B and C resulted in conflicting-but individually strongly supported-relationships among Porphyrio, Fulica, Gallinula, Gallicrex, Amaurornis, Zapornia, Rallina, and Coturnicops (Figures 2 and 3). These conflicting results highlight the need for more complete molecular datasets in order to fully resolve the species-level phylogeny of rails.
Overall, the tree we obtained from our analyses of alignment C closely resembled the results of other phylogenetic studies [27,33,55,[76][77][78][79], although we did identify some important differences. Firstly, while other phylogenetic studies suggested a sister relationship between the clades comprising Amaurolimnas-Aramides-Neocrex-Pardirallus (neotropical wood rails and allies) and Fulica-Gallinula-Porzana-Tribonyx [33,75], we found the Amaurolimnas group to be more closely related to Canirallus (Figure 3). Secondly, we failed to recapitulate two close relationships previously identified by García-R et al. [33]-between Diaphorapteryx hawkinsi (Hawkins's Rail) and Habroptila wallacii (Drummer Rail), and Cabalus modestus (Chatham Islands Rail) and Gallirallus lafresnayanus (New Caledonian Rail)-among the well-studied Pacific rail complex (Gallirallus, Hypotaenidia, and relatives). Thirdly, the results of our analyses strongly supported G. australis (Weka) as sister-taxon to a clade comprising Eulabeornis-Hypotaenidia (Chestnut Rail, and Buff-banded Rail and allies) (Figures 2  and 3), in agreement with some previous studies [55,79], but not with others [33]. Finally, we also found no support for a monophyletic relationship between G. lafresnayanus, G. calayanensis (Calayan Rail), and G. australis (the Gallirallus type-species), a result in agreement with Garcia-R et al. [33] but in conflict with Kirchman [78], and we consequently recommend that G. lafresnayanus and G. calayanensis should be reclassified. The differences between our results and those of previous studies possibly reflect the relative completeness and increased taxon sampling of our alignments (including our new Canirallus sequence) versus previous datasets.
In addition to conflicting topological results, we also observed instances where specimens sequenced in past studies may have been incorrectly identified. For example, two studies have presented mitochondrial data from the south Asian Brown Crake Zapornia akool [77,78]. Our results using Gong et al.'s [77] data supported a close relationship between Zapornia akool, Amaurornis (bush-hens and allies), and Gallicrex cinerea (Watercock), but contradicted previous results based on Ruan et al.'s [78] data, which suggested Z. akool was nested within Zapornia [33]. Morphological data strongly support the former result, and consequently the Z. akool sequences of Ruan et al. are likely to be spurious. As a result, we recommend reinstating the original classification of Z. akool as Amaurornis akool (see Dickinson et al. [80]). Similarly, sequences of Megacrex inepta from the studies of Trewick [81] and those of Kirchman [79] alternatively form incompatible sister-taxon relationships with either G. lafresnayanus (a very close relationship, 98-99% pairwise identity) or the Nkulengu Rail Himantornis haematopus (a relatively distant relationship), respectively ( Figure 3). Morphometric data better support a relationship between H. haematopus and M. inepta [22], suggesting that the sequences from Trewick [81] may originate from G. lafresnayanus. Finally, we also observed that the two Hypotaenidia rovianae sequences included in our analysis were not each other's closest relatives. However, this result may more likely reflect hybridisation, incomplete lineage sorting, or a need for taxonomic revision rather than specimen misidentification (Figure 3).

Morphological Character Evolution
Our reanalysis of Livezey and Zusi's (2007) morphological character matrix (456 characters) revealed three unique (CI = 1.0) and 34 additional non-unique (CI < 1.0) osteological apomorphies uniting Aptornis with Gruiformes. Twelve unique and 22 additional non-unique osteological apomorphies united Aptornis with Heliornis (the only other heliornithoid in the dataset). Our reanalysis of Livezey's [22] morphological character matrix (470 characters) revealed that Aptornithidae was differentiated from Sarothrura-Mentocrex by 36 unambiguous apomorphies. However, there were only 11 apomorphies uniting Aptornis and Sarothrura-Mentocrex, all of which were ambiguous. De Pietri and Mayr [82] found that the Miocene ralloid Paraortygometra porzanoides bears a resemblance to Sarothrura in the morphology of the humerus and tarsometatarsus and suggested a sister relationship between the Heliornithidae and Sarothruridae. We therefore suggest that the characters of the humerus and tarsometatarsus identified by De Pietri and Mayr [81] could be investigated in adzebills as potential unambiguous apomorphies of this newly recognised clade, though as adzebill humeri were vestigial the latter may be more informative.

Node Age Estimates
Node age estimates from our molecular dating analysis are summarised in Table S6. In general, our estimates were similar to or slightly younger than those found by other studies focusing on gruiform birds [27,33,55,56]. However, our mean divergence date estimates between orders and families were slightly older than those found by several recent large-scale phylogenetic studies of birds using genomic data [24,29,30]. These differences likely reflect choice of fossil calibrations, such as use of the Late Cretaceous fossil Vegavis iaai to constrain the divergences of Galloanseres [24,29] or limited use of internal gruiform calibrations [30]. Our older dates may also reflect our exclusive use of mitochondrial loci, which may produce older dates for deeper nodes than genomic data [83]. Nonetheless, our Highest Posterior Densities (HPDs) overlapped considerably with those of Claramunt and Cracraft [24], whose empirically determined priors we used to constrain the age of crown Gruiformes and the split between Procellariiformes and Sphenisciformes.
The age of the split between Aptornis and Sarothrura-Mentocrex was estimated to have occurred 39.6 Ma (25.5-53.8 Ma 95% HPD), during the Middle Eocene. The clade comprising Aptornis and Sarothrura in turn appears to have diverged from the Heliornithidae around 47.8 Ma (34.35-62.2 Ma 95% HPD), during the Middle Eocene. Compared to the heliornithoid clade, extant Rallidae appeared to have a slightly younger crown age of 40.5 Ma (28.4-53.9 Ma 95% HPD). The most recent common ancestor of all Ralloidea (Rallidae and heliornithoids) was estimated to have occurred 54.7 Ma (41.3-70 Ma 95% HPD), which means many fossil Ralloidea (such as messelornithids) [84][85][86], may be young enough to fall within crown Ralloidea. The crown age estimates for extant Gruiformes at 65.3 Ma (49.9-80.4 Ma 95% HPD) and Gruoidea at 59.1 Ma (43.6-76.9 Ma 95% HPD) occurred during the Paleocene to Early Eocene. Finally, within Gruidae we estimated that the divergence between Balearicinae and Gruinae occurred around 24.8 Ma (11.9-39.8 Ma 95% HPD), during the Late Oligocene, and the divergence between Grus and Leucogeranus was estimated at 10.

Discussion
Our phylogenetic analyses unambiguously demonstrated that adzebills (Aptornis; Aptornithidae) are crown gruiforms (Figures 1-3), and we reject the hypothesis that adzebills are close relatives of the Kagu (Rhynochetos jubatus; Eurypygiformes) [21][22][23][24]. Instead, adzebills were found to form an early branch (~39.6 Ma) within Ralloidea as sister-taxon to the Sarothruridae, which we suggest comprises both the flufftails (Sarothrura) and the Madagascan wood rails (Mentocrex). The Aptornithidae-Sarothruridae clade was suggested by our analysis to form a clade with the finfoots (Heliornithidae), comprising the sister-group of the "true" rails (Rallidae) (Figures 1-3). This unexpectedly close relationship between some of the largest and smallest ralloids-adzebills and flufftails, respectively-is not the only instance of discordance between morphological and molecular based phylogenies of Ralloidea [22,33,55,75,88]. For example, Sarothrura was long considered a member of Rallidae, and Sarothrura ayresi was originally described in the Rallidae genus Coturnicops [22,88,89]. In contrast, our analyses confirm previous suggestions that the African Grey-throated Rail (Canirallus oculeus) and the Madagascan wood rails (Mentocrex kioloides sspp.) are not closely related [86,89,90], with the former in fact being a divergent member of Rallidae while the latter was closely related to Sarothrura. The presence of species with "rail-like" morphology in both Rallidae and among heliornithoids (the name we give to the Aptornithidae-Sarothruridae-Heliornithidae clade) suggests that this represents the ancestral state for ralloids, including the ancestors of adzebills.
If the ancestors of heliornithoids resembled archetypal rails, this makes the morphological evolution of adzebills and flufftails, as well as the related finfoots (which are highly-specialised for foot-propelled diving) [91], all the more remarkable. It is also notable that heliornithoids are arguably more morphologically and ecologically varied than the related Rallidae, despite a much lower species diversity [22,90]. All flufftails are diminutive, weighing only~25-50 g [89], whereas Madagascan wood rails may weigh up to~280 g [92], a more typical size for a rail. By comparison, the three species of finfoot range between 120-879 g [92]. In contrast, adult adzebills are estimated to have weighed 16-19 kg (depending on the species), perhaps reaching up to 25 kg [2,3]. The fossil adzebill A. proasciarostratus from the Miocene St Bathans deposits of New Zealand (16)(17)(18)(19) Ma, the only major terrestrial Cenozoic fossil site in New Zealand outside of the Late Quaternary) was only slightly smaller than the two recent species [14]. Thus, if Madagascan wood rails are used as a proxy for the size of the ancestor of adzebills, this suggests that the mass of adzebills increased >50-fold in~20-24 Ma (the time between their estimated common ancestor with Sarothruridae-39.6 Ma-and the age of the St Bathans Fauna). If flighted ancestors of adzebills only arrived in New Zealand following the marine inundation of New Zealand during the Oligocene (peaking during the Waitakian 22-25 Ma) [93] then the temporal window for their size increase would be even further compressed. Similarly rapid evolution appears to have resulted in the presence of flightless rails on numerous oceanic islands (especially taxa in the Gallirallus-Hypotaenidia complex, or in the genera Fulica, Gallinula, Gallirallus, and Porphyrio), many of which became quite large (e.g. the largest Quaternary Rallidae species-NZ's extant Porphyrio hochstetteri and extinct P. mantelli-exceed or exceeded 3-4 kg) [76,81,94,95]. However, none of these rails rival the adzebills for size suggesting either that the adzebill ancestors may have exploited some unique circumstances or that insufficient time has elapsed for extant flightless taxa to reach comparable size (as many extant flightless lineages originated only during the Pliocene or Pleistocene) [76,81,94,95].
We had difficulty identifying morphological features that reliably distinguished heliornithoids from 'true' rails, which is perhaps unsurprising if the more derived lineages-such as the enormous adzebills and the aquatic finfoots-evolved independently from rail-like ancestors. Further, using morphological data to establish the sister taxon of adzebills is problematic due to the antiquity of their flightlessness and their unique feeding specialization, resulting in an "extreme autapomorphy" [22]. By enforcing the Aptornis-Mentocrex-Sarothrura clade in our reanalysis of Livezey's [22] morphological dataset, we showed that no extant rails were recovered as likely "heliornithoids" (including Canirallus). However, it is plausible that other "rail" taxa may yet be revealed as heliornithoids rather than rallids. Genera currently lacking molecular data include Gymnocrex (East Indonesia and New Guinea) and Rougetius (East Africa), both of which have been considered early-diverging members of the Rallidae based on phylogenetic analyses of morphological data [22,88,90]. In addition, Rallicula from New Guinea, which also lacks molecular data, has been simultaneously considered a close relative of Sarothrura and Rallina (incompatible with molecular phylogenies that suggest Sarothrura and Rallina are only distantly related). All other unsequenced extant or recently extinct rail genera are island forms which are likely to have recent dispersal origins (e.g. Aphanapteryx, Cyanolimnas, Mundia, and Pareudiastes) [32]. Our results suggest that substantial work remains to refine the phylogeny of ralloids and that the automatic assignment of species to Rallidae based on "rail-like" morphology is questionable. The latter unfortunately means it may be impossible to correctly discern fossil rails from fossil heliornithoids, which complicates efforts to determine the early biogeography of Ralloidea and especially the geographical source of the adzebill lineage.
Our study reveals a major discordance between the phylogeny and biogeography of the adzebills and their relatives. Unless unrecognised representatives of Sarothruridae occur among extant "rails", the closest living relatives of adzebills are exclusively Afro-Madagascan. Interestingly, this disjunct relationship echoes the sister-taxon relationship observed between the kiwi (Apterygiformes) and Madagascan elephant birds (Aepyornithiformes) [41] and also between the Madagascan Teal (Anas bernieri) and the New Zealand teals (A. aucklandica, A. chlorotis, A. nesiotis, and A. chathamica) [42].
Like these examples, the relationship between the adzebills and Sarothruridae is unlikely to result from continental vicariance, since both the age of the divergence between Aptornithidae and Sarothruridae and the divergence between their common ancestor and Heliornithidae are much too recent (25.5-53.8 and 34.35-62.2 Ma, respectively), versus the separation of Africa and Madagascar from Gondwana >100 Ma [96]. Instead, adzebills likely descended from an ancestral heliornithoid that dispersed to New Zealand overwater from another Gondwanan landmass. It is possible that this ancient dispersal occurred via Antarctica, as geological and palaeontological evidence suggest that at least some coastal regions of Antarctica (and nearby offshore islands) were unglaciated and experienced a temperate climate that supported southern beech (Nothofagaceae) forests until the end of the Eocene (reviewed by Askin and Spicer [95]) or perhaps even the Early Oligocene (see Cantrill [97]). The ancestors of kiwi and moa may also have arrived in New Zealand via Antarctica, as molecular dating results suggest they diverged from their respective nearest living relatives during the Eocene (e.g. Mitchell et al. [41]). In any case, it is likely that "heliornithoids" were formerly more widespread and have subsequently become extinct across much of their former range (with the possible exception of the finfoot lineage, which presently has a pantropical distribution). For example, the Miocene ralloid Paraortygometra porzanoides from France was suggested to have close affinities to Sarothrura [82]. It is possible that the highly derived morphology (and concomitant specialised evolutionary niche) of the adzebills may have contributed to its long persistence, preserving it from whatever environmental conditions promoted the high turnover of more "typical" ralloid lineages observed in the NZ fossil record [98].
We estimated that the North Island and South Island adzebills only diverged relatively recently (0.2-2.3 Ma, with a mean estimate of 1.1 Ma). As the geological precursors of New Zealand's North and South Islands were separated by the Manawatu Strait between~30 Ma and~2 Ma [99], the age of the North and South Island adzebill lineages is incompatible with long-term endemism on both islands. Conversely, the divergence between the two recent adzebill species coincides closely with the formation of an isthmus (1.5-2 Ma) across the Manawatu Strait, which persisted until the modern Cook Strait formed around~0.45 Ma [99,100]. Thus, it is likely that adzebills survived in the larger South Island of New Zealand-from which the fossil adzebill A. proasciarostratus has been described (16)(17)(18)(19) Ma)-and only dispersed into the North Island after the Manawatu land-bridge formed. A similar model has been suggested for several moa species, with the Pleistocene opening and closing of land connection between the North and South Island driving the divergence (1.45 ± 0.8 Ma) between North Island Giant Moa and South Island Giant Moa (Dinornis novaezealandiae and D. robustus, respectively), and between the North Island Mantell's Moa and South Island Heavy-footed Moa (Pachyornis geranoides and P. elephantopus, respectively). As for moa, it is possible that any adzebill lineages endemic to the North Island during the Oligocene subsequently became extinct during the marine transgression (peaking around 23 Ma) that heavily reduced the land area of what would eventually become the North Island [99]. South Island adzebills were slightly larger on average than North Island adzebills, a pattern common in other New Zealand birds with North/South Island sister-species (such as geese in Cnemiornis or moa in Pachyornis) and is consistent with Bergmann's rule [22]. Otherwise, both species were physically near-identical. It is therefore likely that adzebills occupied the same niche in the North Island as they did in the South Island, and apart from a slight change in mass, little to no macroevolution occurred since the species diverged.
In conclusion, our molecular data revealed adzebills to be crown ralloids, contrary to a long-standing hypothesis that they were close relatives of the New Caledonian kagu. The resolution of this centuries-old taxonomic issue implies that adzebills were not of Gondwanan vicariant origin, but instead descended from a rail-like bird that arrived in New Zealand by long-distance overwater dispersal during the latest Eocene. Consequently, the ancestors of the adzebill almost certainly became flightless after their arrival in New Zealand. Our new results-confidently placing the giant adzebills within ralloids-further reinforce Ralloidea as an ideal clade for studying the genomic drivers and consequences of flightlessness. The independent loss of flight in adzebills over 16 million years ago stands in stark contrast to the more recent (Pliocene/Pleistocene) losses of flight in extant island rail species.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1424-2818/11/2/24/s1, Figure S1: MapDamage report for the final round of mapping reads from DM8046 against the DM8046 majority-rule consensus. The top panels show the characteristic high frequency of purines (A and G) immediately before the reads. The two lower panels show the accumulation of 5' C-to-T (red) and 3' G-to-A (blue) misincorporations characteristic of ancient DNA; Figure S2: MapDamage report for the final round of mapping reads from S23033 against the DM8046 majority-rule consensus. The top panels show the characteristic high frequency of purines (A and G) immediately before the reads. The two lower panels show the accumulation of 5' C-to-T (red) and 3' G-to-A (blue) misincorporations characteristic of ancient DNA. Figure S3: MapDamage report for the final round of mapping reads from CM AV18536 against the DM8046 majority-rule consensus. The top panels show the characteristic high frequency of purines (A and G) immediately before the reads. The two lower panels show the accumulation of 5' C-to-T (red) and 3' G-to-A (blue) misincorporations characteristic of ancient DNA; Figure S4: Saturation plots of the primary alignment (predominantly complete mitochondrial genomes) generated using DAMBE. Blue points represent transitions (s), green points represent transversions (v).; Table  S1: Details of new mitochondrial genome sequences; Table S2: Mitochondrial genome sequences used for bait design for enrichment of the CM AV18536 eggshell library; Table S3: List of all sequences used in this study with GenBank accession numbers; Table S4: Saturation test statistics as generated in the test by Xia et al. (complete mitochondrion genome sequences only); Table S5: Partitioning schemes used for all phylogenetic and molecular clock analyses; Table S6: Age estimates for key nodes from our BEAST analyses.