Erschienen in: Molecular Phylogenetics and Evolution ; 101 (2016). - S. 303-313
https://dx.doi.org/10.1016/j.ympev.2016.05.021

A tree of geese: A phylogenomic perspective on the evolutionary historyof True Geese
Jente Ottenburghs a,⇑, Hendrik-Jan Megens b, Robert H.S. Kraus c,d, Ole Madsen b, Pim van Hooft a,
Sipke E. van Wieren a, Richard P.M.A. Crooijmans b, Ronald C. Ydenberg a,e, Martien A.M. Groenen b,
Herbert H.T. Prins aaResource Ecology Group, Wageningen University, Droevendaalsesteeg 3a, 6708 PB Wageningen, The Netherlands
Animal Breeding and Genomics Centre, Wageningen University, Droevendaalsesteeg 1a, 6708 PB Wageningen, The NetherlandscDepartment of Biology, University of Konstanz, D-78457 Konstanz, GermanydDepartment of Migration and Immuno-Ecology, Max Planck Institute for Ornithology, Am Obstberg 1, D-78315 Radolfzell, GermanyeCentre of Wildlife Ecology, Simon Fraser University, V5A 1S6 Burnaby, BC, Canadaba b s t r a c t

Keywords:
Consensus
Concatenation
Gene tree
Hybridization
Incomplete lineage sortingSpecies tree

Phylogenetic incongruence can be caused by analytical shortcomings or can be the result of biologicalprocesses, such as hybridization, incomplete lineage sorting and gene duplication. Differentiation
between these causes of incongruence is essential to unravel complex speciation and diversiﬁcation
events. The phylogeny of the True Geese (tribe Anserini, Anatidae, Anseriformes) was, until now, contentious, i.e., the phylogenetic relationships and the timing of divergence between the different goose
species could not be fully resolved. We sequenced nineteen goose genomes (representing seventeen spe
cies of which three subspecies of the Brent Goose, Branta bernicla) and used an exon based phylogenomic
approach (41,736 exons, representing 5887 genes) to unravel the evolutionary history of this bird group.
We thereby provide general guidance on the combination of whole genome evolutionary analyses andanalytical tools for such cases where previous attempts to resolve the phylogenetic history of several taxa
could not be unravelled. Identical topologies were obtained using either a concatenation (based upon analignment of 6,630,626 base pairs) or a coalescent based consensus method. Two major lineages, corre
sponding to the genera Anser and Branta, were strongly supported. Within the Branta lineage, the White
cheeked Geese form a well supported sub lineage that is sister to the Red breasted Goose (Branta ruﬁcollis). In addition, two main clades of Anser species could be identiﬁed, the White Geese and the Grey Geese.
The results from the consensus method suggest that the diversiﬁcation of the genus Anser is heavily inﬂu
enced by rapid speciation and by hybridization, which may explain the failure of previous studies toresolve the phylogenetic relationships within this genus. The majority of speciation events took placein the late Pliocene and early Pleistocene (between 4 and 2 million years ago), conceivably driven by aglobal cooling trend that led to the establishment of a circumpolar tundra belt and the emergence of tem
perate grasslands. Our approach will be a fruitful strategy for resolving many other complex evolutionary
histories at the level of genera, species, and subspecies.

1. Introduction
Incongruence between phylogenies generated from different
sets of genetic and phenotypic data poses a signiﬁcant challenge

Abbreviations: BDS, birth–death process with species sampling; bp, base pairs;
CI, conﬁdence interval; cytb, cytochrome b; GTR, general time reversible; mtDNA,
mitochondrial DNA; Mya, million years ago; ND2, NADH dehydrogenase 2.
⇑ Corresponding author.
E-mail address: jente.ottenburghs@wur.nl (J. Ottenburghs).

in evolutionary biology. Such incongruence can be caused by ana
lytical shortcomings (Davalos et al., 2012; Rokas et al., 2003), suchas issues with limited taxon sampling, unmet assumptions in themodelling of sequence evolution, and the choice of different opti
mality criteria (Graybeal, 1998; Rokas et al., 2003; Yang et al.,1994). Phylogenetic incongruence can also be generated by biolog
ical processes, such as hybridization, incomplete lineage sortingand gene duplication (Degnan and Rosenberg, 2009; Maddison,
1997; Pamilo and Nei, 1988). Differentiation between these causes

Konstanzer Online-Publikations-System (KOPS)
URL: http://nbn-resolving.de/urn:nbn:de:bsz:352-0-379362

304

of incongruence is essential to unravel complex speciation anddiversiﬁcation events.
Although phylogenetic incongruence compromises the estimation of species trees, it can be a virtue for the evolutionary biolo
gist, who attempts to better understand complex speciation
events that are often obscured by rapid diversiﬁcation andhybridization (Degnan and Rosenberg, 2009; Toews andBrelsford, 2012). Genomes are mosaics of many different gene his
tories (Ellegren et al., 2012; Maddison, 1997; Payseur, 2010) thatcan be compared to deduce the species tree, while also gaining
more insights into intricate events of the speciation process, suchas hybridization episodes or bursts of rapid diversiﬁcation (Innan
and Watanabe, 2006). In addition, conﬂicting gene histories canalso be used to infer population sizes and divergence times(Rannala and Yang, 2003; Wall, 2003). Analysing phylogenetic
incongruence from a genomic perspective of diverse gene histories
will lead to a more complete picture of the speciation process.
The occurrence of phylogenetic incongruence is a common andwidespread phenomenon: analyses of different genes often lead to
discordant gene trees (Rokas et al., 2003). Consequently, several
methods have been developed to estimate the species tree froma collection of contrasting gene trees (Degnan and Rosenberg,
2009; Delsuc et al., 2005). In concatenation (or supermatrix) methods all sampled genes are concatenated and analysed as a single‘‘supergene” (de Queiroz and Gatesy, 2007). Supertree methods,
on the other hand, involve separate analyses of the sampled genes
and subsequent integration of the resulting trees into a species tree
(Bryant, 2003). Certain supertree methods incorporate the multi
species coalescent model to estimate the species tree from a setof heterogeneous gene trees (Knowles, 2009). The multispecies
coalescent model extends the classical coalescent (Kingman,1980) to multiple populations and describes gene trees as indepen
dent random variables generated from the coalescence process
occurring along lineages of the species tree (Liu et al., 2015). A
number of coalescent based methods have been developed to esti
mate the species tree from multigene sequence data, including
Bayesian (e.g., BUCKy, Ane et al., 2007; BEAST, Drummond et al.,2012) and likelihood approaches (e.g., STEM, Kubatko et al.,2009; STELLS, Wu, 2012). Coalescent based methods have beensuccessfully implemented to disentangle the complex evolutionary
history of several closely related bird species (Carling et al., 2010;
Hung et al., 2012; Lavretsky et al., 2014).
The failure to unravel the phylogenetic relationships between
goose species of the tribe Anserini (commonly referred to as True
Geese, see Table 1 for a taxonomic overview of this tribe) can be
attributed to high levels of phylogenetic incongruence (Ruokonen
et al., 2000). Although there have been several studies delving into
the phylogeography of particular goose species (Humphries et al.,2009; Jonker et al., 2013; Paxinos et al., 2002; Quinn, 1992;
Ruokonen et al., 2005, 2008; Scribner et al., 2003; Volkovsky
et al., 2013), the phylogenetic relationships between most species
remained unclear. Until now, species level phylogenies were based
either on morphology (Livezey, 1996), or mitochondrial markers(Donne Gousse et al., 2002; Gonzalez et al., 2009; Lee et al.,2008; Ruokonen et al., 2000). These approaches resulted in contra
dicting topologies and different mitochondrial markers led to dif
ferent phylogenetic reconstructions (Fig. 1). Phylogenetic
incongruence is particularly apparent in the genus Anser. Although
phylogenetic incongruence can be caused by hybridization and/or
incomplete lineage sorting, all previous studies on goose phylogeny were based on mitochondrial DNA (mtDNA), which sug
gests analytical issues as the most likely cause of theseincongruences. One example of the discrepancies among previous
studies involves the phylogenetic position of the Lesser White
fronted Goose (Anser erythropus). Morphological analyses(Livezey, 1996) and some mitochondrial markers (cytb and ND2,

Table 1
Current taxonomy for the True Geese (tribe Anserini) based on Clements et al. (2015).
English nameScientiﬁc name

Genus ANSER
Swan Goose
Taiga Bean Goose

Anser cygnoidesa
Anser fabalis

Tundra Bean Goose

Anser serrirostris

Pink-footed Goose

Anser
brachyrhynchusa
Anser albifrons

Greater White-fronted
Goose

A. f. fabalisa
A. f. johanseni
A. f. middendorfﬁi
A. serrirostris rossicusa
A. serrirostris serrirostris

A. a. albifrons (Eurasian)aA. a. ﬂavirostris (Greenland)
A. a. gambeli (Western)
A. a. frontalis (Western)
A. a. elgasi (Tule)

Anser erythropusa

Lesser White-fronted
Goose
Greylag Goose

Anser anser

Bar-headed Goose
Emperor Goose
Snow Goose

Anser indicusa
Anser canagicusa
Anser caerulescens

Ross’ Goose

Anser rossiia

Genus BRANTA
Brent Goose

Branta bernicla

Barnacle Goose
Cackling Goose

Branta leucopsisa
Branta hutchinsii

Canada Goose

Branta canadensis

Hawaiian Goose

Brantasandvicensisa
Branta ruﬁcollisa

Red-breasted Goosea

Subspecies

A. a. anser (European)aA. a. rubrirostris (Siberian)A. c. caerulescensa
A. c. atlantica

B. b. bernicla (Dark-bellied)aB. b. hrota (Pale-bellied or
Atlantic)aB. b. nigricans (Black)aB. b. orientalis
B. h. leucopareia (Aleutian)B. h. hutchinsii (Richardson’s)
B. h. minima (Minima)aB. h. taverni (Taverner’s)
B. c. mofﬁtti
B. c. maxima
B. c. occidentalis
B. c. fulva
B. c. canadensisa
B. c. interior
B. c. parvipes

Blood samples retrieved from this (sub)species.

Donne Gousse et al., 2002; Gonzalez et al., 2009) report a sister
relationship with the Greater White fronted Goose (Anser alb
ifrons), whereas studies using the mitochondrial control region(Lee et al., 2008; Ruokonen et al., 2000) place the Lesser White
fronted Goose in a clade with Taiga Bean Goose (Anser fabalis)and Pink footed Goose (Anser brachyrhynchus).
In accordance with the lack of a well supported phylogeny, thetiming of origin of the extant goose species is unknown too. Fossil
evidence indicates that geese were present during the Miocene andPliocene (Brodkorb, 1964) and several phylogeographic studies
reported Pleistocene origins of certain goose subspecies (Avise
et al., 1992; Shields, 1990; Van Wagner and Baker, 1990). Moreover, a mtDNA study of the genus Anser dated speciation eventsto the late Pliocene and early Pleistocene (Ruokonen et al., 2000).
Based on the available evidence, we hypothesize that the diversiﬁ
cation of modern goose species was initiated in the Pliocene,
accompanied by further diversiﬁcation into distinct subspecies
during the Pleistocene. Recent developments in divergence timing
methods, that integrate both genomic data and fossil calibration
points, can be applied to test this hypothesis (Parham et al., 2012).
The objectives of this study were threefold: (1) to unravel thephylogenetic relationships within the Anserini tribe using phyloge
nomic tools for species tree estimation, (2) to assess the timing of

305

A. albifrons
A. erythropus
A. anserA . fabalis
A. brachy rhynchus
A. rossii
A. caerulescens
A. canagicus
A. lndicus
A. cygnoides
B. c anadensis
B. sandvicensls
B. leucopsls
B. bemicla
B. ruficolis

Livezey (1996) Morpho logy

B

A. f abais
A. brachyrhynchus

A . al:>~rons

A. erythropus
A. abfrons

A . f abais
A . anser

A. anser
A. rossi

A . cygnoides
A . brachyrhynchus

A. lndicus

A .lndicus

A. canagicus
A. caerulescens

A . rossi
A . canagicus

B. canadensis
B. leucops Is
B. ruficolis

B. canadensis
B. leucopsis
B. rufocolis

B. bernicla

B. bernicla

Lee et al. (2008) mtDNA control region

A . erythropus

Gonzalez et al. (2009) cytb & ND2

Fig. 1. Phylogenetic trees for the True Geese based on morphology or mitochondrial markers. (A) Morphology ( Livezey. 1996). ( B) Mitochondrial control region ( Lee eta!..
2008}. ( C) Cytochrome b and NAOH hydrogenase 2 ( Gonzalez et aL. 2009}

dive~gence for the extant goose species by means of genomic andfossil data and (3) discuss these findings in a framework of ecolog
ical, biogeographic and climatic events.

22. Processing sequences(representing 5887 genes) were concatenated, resulting in analignment of 6,630,626 base pairs (bp ).
Paired end reads {1 00 bp) were also mapped to the Swan Goose(Anser cygnoides) genome (Lu et a l, 201 5) using BWA ( Li andDurbin, 2009). Over 95% of the reads were mapped successfully
for all species. However, mapping to the genome of a species thatis included in the phylogenetic analysis (in this case the Swan
Goose) potentially leads to an inherent bias that could interferewith the phylogenetic analysis. For instance, mean mapping qual
ity (Fig. 52) is higher for species belonging to the same genus(Anser) as the reference genome compared to species of the othergenus (Bmnta). Therefore we decided to use the variant calls basedon the mapped reads to a more distantly related reference genome,
the Mallard (A platyrhynchos). Although this choice leads to a
reduction in the amount of data, it avoids the bias introduced by
mapping to the Swan Goose genome, because all species areequally distantly related to the Mallard. A similar approach has
been used to reconstruct the phylogeny of Cichlid fish {lives and
Lopez Fernandez, 2014).

Paired end reads {lOObp) were mapped to the Mallard (Anasplatyrhynchos) genome, version 73 (Huang et al., 2013, supplemen

2.3. Phylogenomic analyses

tary material) using SMALT (http://www.sanger.ac.uk/resourcesfsoftwarefsmalt/). Previous research indicated that SMALT
is appropriate for mapping paired end reads when the reference
genome is distantly related to the sampled species (Frantz et al.,2013 ). Furthermore, preliminary analyses showed that this soft
ware produces the best results compared to other mapping soft
ware (Fig. 52). SMALT uses a hash index of short words, sampled
at equidistant steps along the reference genome. We applieddefault settings (word length; 13, spacing of hashed words; 6 ).
Over 99% of the reads were mapped successfully in all samples,but to decrease the incidence of off site mapping we only accepted
properly mapped paired reads, resulting in mapping rates between
63% and 78% (Table 52 ~ From the resulting bam files duplicate
sequences were excluded using SAMtools dedup and realignedwith lndeiRealigner in GATK 2.6. Variant sites were called bymeans of the UnifiedGenotyper in GATK 2.6 with a minimum basequality of20 and a heterozygosity value ofO.Ol (this value has pro
ven to be suitable for variant calling after mapping to a distantly
related reference genome, see Frantz, 2015 ). The genomic positionsof exonic sequences that were one to one orthologous between
the Mallard and other bird species (chicken, turkey, flycatcher
and zebra finch ) were retrieved from the ENSEMBL data base.
These genomic positions were used to extract the sequences fromthe different goose genomes. A final set of 41,736 unique exons

We applied several phylogenomic methods based on maximum
likelihood implemented in RAxML 8.3 (Stamatakis, 2006 ). We usedboth supertree and supermatrix techniques (Delsuc et al., 2005).
First, the concatenated alignment was analysed under two com
monly used substitution models (GTR + r and GTR + r + 1), bothwith 100 bootstrap replications (Ranwez et at., 2007). Second, weselected 3570 one to one orthologous genes {thus combining thegene exons) with a minimum length of 500 bp. These genes wereanalysed separately under a GTR + r substitution model with
100 rapid bootstraps. The resulting gene trees were filtered on
average bootstrap support (minimum >50). This final set of 3558
well supported gene trees was consequently combined into a con
sensus tree using the software package STELLS version 1.6.1 (Wu,2012). STELLS infers the species tree from a given set of gene trees,using the multispecies coalescent as underlying genealogical pro
cess. To assess support for particular clades in the supertree analysis, concordance factors were calculated in DensiTree ( Bouckaert,
2010). A concordance factor represents the percentage of genetrees that contain the same node as the species tree (Baum,2007; Knowles and Kubatko, 2010). Trees were rooted either byusing an outgroup (Mallard), or by midpoint rooting (without anoutgroup). Both methods resulted in the same position of the rootin the tree. Here, we present the results of the midpoint rooting

2. Materials and methods
2.1. Sample collection and DNA extraction

We collected blood samples from 19 goose (sub )species (Table 1,
see Table 51 for details). Genomic DNA was isolated from the blood
samples using the Qiagen Gentra kit (Qiagen Inc.). Quality andquantity of the DNA was measured using the Qubit {Invitrogen, Life
Technologies). Sequence libraries were made according to lllumina
protocols and sequenced paired end {100 bp) on the Hiseq2500(lllurnina Inc.).

306

only because the high divergence of the root species (Mallard)
resulted in conﬂicting parameter estimations in RAxML 8.3.

2.4. Molecular clock analyses
Ideally, a molecular clock analysis is run with multiple calibration points (fossils or biogeographic events) using a relaxed clock(Parham et al., 2012). Unfortunately, the fossil record for geese is
too sparse to use multiple calibration points (see Section 4). Therefore, we followed the method outlined below. The resulting divergence times should thus not be interpreted as the exact dates ofdivergence, but rather as rough guiding estimates that allow us
to formulate testable hypotheses on the biogeographical and eco
logical drivers of goose speciation.
We estimated divergence times using an approximate likeli
hood method as implemented in MCMCtree (in PAML version 4,
Yang, 2007), with a global clock and birth death sampling. We ﬁt
ted a GTR + C model to the concatenated alignment and estimateda mean substitution rate by ﬁtting a strict clock and setting theroot age at 9 Mya based on previous studies (Fulton et al., 2012;
Jetz et al., 2012). This mean substitution rate (0.134992 ±0.000307 substitutions/site/Myr) was used to adjust the prior on
the substitution rate (rgene) modelled by a C distribution as C
(1, 7.4). Parameters for the birth death process with species sampling (BDS) and r2 values were set at 1 1 0 and G(1, 10), respec
tively. The Anser Branta divergence was constrained between 20
and 4 million years ago to reﬂect reports of goose fossils in thistime period (Brodkorb, 1964). MCMCtree analyses were run multi
ple times to check for convergence of results. Analyses with arelaxed clock gave similar results compared to the global clock
approach, but without convergence of results. Therefore, we only
present the estimates of the global clock analysis here. To further
assess the reliability of the MCMCtree analysis, we performed a
global clock calculation of divergence times using the software
package PATHd8 (Britton et al., 2007), setting the root age at
9.5 Mya as a result of the MCMCtree analysis.

3. Results
The concatenation and consensus methods resulted in identical
topologies (Fig. 2). All clades in the concatenation analysis were
supported by maximum bootstrap values, except for the sister
species relation between Tundra Bean Goose and Pink footed
Goose which has a bootstrap support of 95. Analyses with the
GTR + C substitution model and the GTR + C + I substitution model
resulted in the same topology. There is a clear separation between
the two genera (Anser and Branta). The genus Branta consists of aclade of White cheeked Geese Canada Goose (Branta canadensis),
Cackling Goose (Branta hutchinsii), Barnacle Goose (Branta leucopsis) and Hawaiian Goose (Branta sandvicensis) and two basal splits
leading to Brent Goose (Branta bernicla) and Red breasted Goose(Branta ruﬁcollis). The most basal split in the genus Anser leads to
the Bar headed Goose (Anser indicus). Next, two main clades canbe recognised: the White Geese Snow Goose (Anser caerulescens),
Ross’ Goose (Anser rossii) and Emperor Goose (Anser canagicus)
and the Grey Geese
Greylag Goose (Anser anser), Swan Goose(A. cygnoides), the White fronted Geese (A. albifrons and A. erythropus) and the Bean Goose complex (A. fabalis, Anser serrirostris and A.
brachyrhynchus). Concordance factors ranged between 4% and 93%(Fig. 2), with generally higher values found within Branta (30 93%)
than within Anser (4 54%). Indeed, concordance factors in the con
sensus analysis were consistently higher for the Branta comparedto the Anser (Mann Whitney U = 60.5, p = 0.015), which indicateshigher levels of phylogenetic incongruence in the latter clade.

Fig. 3 gives an overview of the divergence times (with 95% CI)
for the different speciation events. The split between Anser and
Branta was estimated at 9.5 Mya (15.1 4.2). This estimate wasused for root calibration in the PATHd8 analysis, which resulted
in similar divergence times compared to the MCMCtree analysis(R2 = 0.973, p < 0.001, Fig. S3). Most mean divergence times areconcentrated in between 4 and 2 Mya. The mean divergence times
for Red breasted Goose (B. ruﬁcollis) and Brent Goose (B. bernicla)
are dated to approximately 6.8 and 5.8 Mya, respectively. Finally,
three Brent Goose subspecies (bernicla, hrota and nigricans) origi
nated during the Pleistocene, around one Mya. An overview ofthe evolutionary history of True Geese is presented in Fig. 4, com
bining phylogenetic relationships, divergence times and current
distributions.
4. Discussion
Whole genome analyses have recently been used to unravel thephylogenetic relations between the major bird orders (Jarvis et al.,2014; Prum et al., 2015), but our study presents one of the ﬁrst
phylogenomic analyses of a group of closely related bird species(DaCosta and Sorenson, 2016; Nater et al., 2015). As both concate
nation and consensus methods resulted in the same topology andall nodes on the concatenation tree are supported by high bootstrap values, it is likely that the resulting phylogeny closely
approaches the actual species tree for the Anserini. We ﬁrmly
resolve several incongruences among previous studies (seeSection 4.1).
The resolution of these incongruences can be attributed to thedeliberately exon based phylogenomic approach in our study(Ilves and Lopez Fernandez, 2014). First, the choice of one to one
orthologous exonic sequences ensures high mapping quality
because of greater sequence conservation of exons compared to
introns (Kerr et al., 2014; Ranwez et al., 2007). Moreover, the useof one to one orthologous sequences avoids the analytical com
plexities introduced by duplicated genes (Maddison, 1997) andthe occurrence of indels and greater sequence divergence thatcan pose substantial problems during multiple sequence alignment(Liu et al., 2010). Second, the genome wide set of thousands ofgenes enabled us to deduce with high statistical power the species
tree despite potentially high levels of phylogenetic incongruencecaused by incomplete lineage sorting and/or hybridization. Indeed,both the concatenation and the consensus method uncovered thesame phylogenetic signal from the genome wide sequence data.
Furthermore, the concordance factors calculated in the
coalescent based consensus method provide important insights
into the effects of incomplete lineage sorting and/or hybridizationin certain clades.
4.1. Phylogenetic relationships
Monophyly of the genera Anser and Branta was already wellestablished (Delacour and Mayr, 1945). However, the timing ofseparation of these genera is still a matter of debate. Fossil evi
dence suggested that these two groups of geese have a common
ancestor between 5 and 4 Mya (Brodkorb, 1964; Wetmore, 1956),
and this date has been used to calibrate the molecular clock inbirds (Shields and Wilson, 1987). The consistency of this molecular
clock has been questioned (van Tuinen and Hedges, 2001; Weirand Schluter, 2008), because subsequent molecular studies
reported older dates for the Anser Branta split ranging from 23
to 9 million years ago (Brown et al., 2008; Fulton et al., 2012;
Gonzalez et al., 2009; Jetz et al., 2012; Jiang et al., 2010; Pereiraand Baker, 2006). We decided to use 9 million years to calibrate
the mutation rate, based on recent estimates (Fulton et al., 2012;

307

100 1

100 1

Dar1<-bellied Brent Goose (B. bemicla bemida)

37.621Pale-bellied Brent Goose (B. bemicla hrota)

93.931

Black Brent Goose (B. bemida nigrlcans)
Canada Goose (B. canadensis)

100
91.09
100

100 1
35.481

74.19

100

Barnacle Goose (B. leucopsis)

100 130.42 1

Cackling Goose (B. hutchinsii)
Hawaiian Goose (B. sandvicensls)

48.92

Red-breasted Goose (B. ruficollis)
Bar-headed Goose (A. indicus)
Greater White-fronted Goose (A. albifrons)

100 1
9.53 I

Lesser White-fronted Goose (A. erythropus)

100
4.89
100 1
9.281

100

91.09

8.35

Tundra Bean Goose (A. senirostris)

95 1
13.581

Pink-footed Goose (A. brachyrhynchus)
Taiga Bean Goose (A. fabalis)

100
13.89

Swan Goose (A. cygnoides)
Greylag Goose (A. anser)

100
21.02

l

100
17.081

100

Snow Goose (A. caerulescens)

54.541

Ross· Goose (A. rossii)
Emperor Goose (A. canagicus)

Fig. 2. Maximum Likelihood Tree for the True Geese based on consensus and concatenation methods. Bootstrap values of the concatenation analysis above branches,
concordance factors of the consensus analysis below.

MIOCENE

PUOCENEPlEISTOCENE
B.. b. ~

I

I

I

I

I
I

I

I

B. ~

a. I.._.;.

B.-

I
I
I

I
I
I

I

I
15

12.5

10

7.5

5

2.5

Fig. 3. limed phylogenetic tree using a global clock approach in MCMCtree. Blue bars indicate 95% confidence intervals.

0

308

Fig. 4. An overview of the evolutionary history of the True Geese combining phylogenetic tree, divergence times and current distributions of all species. Distributions basedon BirdLife’s species range maps. Drawings used with permission of the Handbook of Birds of the World.

Jetz et al., 2012). Moreover, we used a fossil constraint between 20
and 4 million years ago to reﬂect reports of goose fossils in thistime period (Brodkorb, 1964). The resulting wide conﬁdence interval for the Anser Branta split in our molecular clock analysis, rang
ing from 15.1 to 4.2 Mya, indicates that there is still considerable

uncertainty for this estimate. This uncertainty is a consequenceof the lack of proper fossil calibration points for the Anserini tribe.
Although there are numerous goose fossils (Brodkorb, 1964), it is
not possible to conﬁdently determine the phylogenetic positionof these fossils. A thorough phylogenetic analysis of goose fossils,

309

Fig. 4 (continued)

such as the one recently published for the Charadriiformes (Smith,2015), is clearly essential in order to achieve better estimates forthe divergence times.
Phylogenetic relationships within the genus Branta (commonlyreferred to as the Black Geese) could be unambiguously inferred,consisting of a clade of White cheeked Geese Canada Goose (B.
canadensis), Cackling Goose (B. hutchinsii), Barnacle Goose (B. leu
copsis) and Hawaiian Goose (B. sandvicensis) and two basal splits
leading to Brent Goose (B. bernicla) and Red breasted Goose (B.
ruﬁcollis). The exact phylogenetic position of the latter two species

310

remained uncertain up to our study: some studies reported a sister
species relation (Lee et al., 2008; Paxinos et al., 2002), while others
placed either Red breasted Goose (Donne Gousse et al., 2002;
Livezey, 1996) or Brent Goose at the base of the Branta clade(Gonzalez et al., 2009). We ﬁnd strong support for a basal split
leading to the Brent Goose in Branta. Within the White cheeked
Geese clade there is a clear distinction between Canada Goose
and Cackling Goose which were considered large and small
bodied forms of the same species in the past (Owen, 1980). The sig
niﬁcant morphological and molecular divergence between these
species reported in previous studies is therefore conﬁrmed byour genomic analysis (Scribner et al., 2003; Van Wagner andBaker, 1990). In addition, we ﬁnd support for the sister species
relation between Cackling Goose and Barnacle Goose described
previously (Paxinos et al., 2002).
The Hawaiian Goose was the ﬁrst species to split from theWhite cheeked Geese clade, around 3.5 Mya (5.5 1.5). This datesuggests that the ancestors of the Hawaiian Goose arrived on
islands that were present in this time period, such as Kauai, Nihoa
or Oahu, (Macdonald et al., 1983). However, it is also possible thatthe divergence of the Hawaiian Goose occurred on the mainland,followed by colonization of Hawaii at a later stage. Our results con
trast with the ﬁndings of Paxinos et al. (2002), who reconstructed
the evolutionary history of the Hawaiian goose radiation (i.e. themodern Hawaiian Goose and several extinct species) using mtDNA
extracted from fossil remains. They found that the Hawaiian radi
ation is a sister group to the Canada Goose that diverged approxi
mately 890,000 years ago (1.3 Mya 577,000), a date also reportedby Quinn et al. (1991) based on restriction fragments of mtDNA.
The incongruence between our analysis and those of Paxinos
et al. (2002) and Quinn et al. (1991) is probably the result of
incomplete lineage sorting. A closer look at the results of the con
sensus method (which takes into account incomplete lineage sorting) supports this explanation. The clade comprising Canada
Goose, Cackling Goose and Barnacle Goose is recovered in 35.5%of the gene trees, whereas a sister relation between Canada Gooseand Hawaiian Goose is supported by 16% of the gene trees. Gene
tree incongruence can also be the result of hybridization(Maddison, 1997), but this is less likely as the Hawaiian Goose
has been geographically isolated from the Canada Goose. However,
we cannot rule out the possibility that vagrant Canada Geese might
have reached the Hawaiian archipelago in the past (Berger, 1972)or that Canada Geese and Hawaiian Geese have hybridized before
the colonization of Hawaii. Another incongruence between our
results and those of Paxinos et al. (2002) and Quinn et al. (1991)
concerns the timing of divergence of the Hawaiian Goose(3.5 Mya versus 890,000 years ago). This large difference can be
explained by the fossil calibration used by these studies as theyboth calibrated the molecular clock by setting the date of theAnser Branta split at 4.5 Mya, which is later than the date we use
here (9.5 Mya). In addition, the effective population size of mtDNA
is smaller compared to nuclear DNA, which results in a consider
ably younger divergence time (Ballard and Whitlock, 2004). There
fore the Hawaiian Goose most likely diverged earlier than
expected, before the further diversiﬁcation of the white cheeked
goose ancestor that gave rise to the Canada Goose, Cackling Gooseand Barnacle Goose.
In contrast to Branta, the phylogenetic relationships within thegenus Anser were largely unresolved up to this study. The failure toreconstruct the evolutionary history of this bird group can be
attributed to rapid speciation, which results in incomplete lineage
sorting, and hybridization (Ruokonen et al., 2000). Indeed, concor
dance factors in the consensus analysis were lower for the Anser
clade compared to the Branta clade, indicating higher levels ofgene tree discordance in the former clade. By using genomic data
and methods that take into account gene tree discordance, wewere able to construct a well supported phylogeny for the genus
Anser. The most basal split leads to the morphologically divergent
Bar headed Goose (A. indicus). Next, two main clades can be recog
nised: the White Geese Snow Goose (A. caerulescens), Ross’ Goose(A. rossii) and Emperor Goose (A. canagicus) that are sometimes
placed in a separate genus, Chen (Gil and Donsker, 2013), and the
Grey Geese Greylag Goose (A. anser), Swan Goose (A. cygnoides),the White fronted Geese (A. albifrons and A. erythropus) and theBean Goose complex (A. fabalis, A. serrirostris and A.
brachyrhynchus).
The monophyly of the White Geese was already wellestablished from a morphological perspective (Livezey, 1996),
but molecular analyses were unable to recover this clade(Donne Gousse et al., 2002; Gonzalez et al., 2009; Lee et al.,2008; Paxinos et al., 2002). We were able to conﬁrm the mono
phyly of this clade, although gene tree discordance was apparent
in the consensus method as only 17.1% of the gene trees supported
this clade. From a taxonomic point of view the White Geese should
be classiﬁed in the genus Anser, because placing them in the dis
tinct genus Chen results in polyphyly of the genus Anser.
Within the Grey Geese clade, the position of the Swan Goose, asa sister species to the White fronted Geese and the Bean Goose
complex, is quite surprising. Given its morphological disparity,
the Swan Goose was traditionally placed at the base in the Grey
Geese clade (Livezey, 1996) and has even been classiﬁed in a sep
arate genus, Cygnopis (Owen, 1980). Also, the phylogenetic positionof the Lesser White fronted Goose has been established as a sister
species of the Greater White fronted Goose. The placement of theLesser White fronted Goose in a clade with Taiga Bean Gooseand Pink footed Goose by previous studies can be explained byanalytical shortcomings, incomplete lineage sorting and/or
hybridization (Lee et al., 2008; Ruokonen et al., 2000).
The taxonomy and evolutionary relationships of the Bean Goose
complex have been controversial (Delacour, 1951; Mooij andZöckler, 1999; Sangster and Oreel, 1996). The current classiﬁcationencompasses three species: the Pink footed Goose, the Taiga Bean
Goose (three subspecies) and the Tundra Bean Goose (two sub
species). However, based on analyses of the mitochondrial control
region, Ruokonen et al. (2008) identiﬁed three distinct lineages:
the Pink footed Goose, the Middendorff’s Goose (currently a sub
species of Taiga Bean Goose), and the Bean Goose (currently split
in Taiga and Tundra Bean Goose). We report a sister species relation between Pink footed Goose and Tundra Bean Goose. The phylogenetic relationships in this complex are highly inﬂuenced byincomplete lineage sorting and/or hybridization. A thorough geno
mic analysis in a phylogeographic context with broad sampling
across the entire range of this species group is needed to conﬁ
dently reconstruct the evolutionary history of the Bean Goose
complex.
4.2. Biogeography and ecology
It can be conceived that True Geese originated in Eurasia, considering the fact that the earliest (Miocene) goose fossils have been
excavated in Western Europe (Brodkorb, 1964) and species at thephylogenetically most basal splits (Bar headed Goose, Brent Gooseand Red breasted Goose) currently have a largely Eurasian distri
bution (although the current distribution does not necessarily
coincide with the past distribution). If we assume a Eurasian origin,
at least three colonization events to North America can be recog
nized: by the White cheeked Geese, by the White Geese, andﬁnally by the Greater White fronted Geese. However, this scenario
remains to be tested by means of ancestral reconstructions in com
bination with fossil data. Based on the phylogeny and present day
biogeographic patterns, we can recognize three main groups: theBlack Geese (Branta) with members in North America and Eurasia,

311

the Grey Geese (Anser) with a largely Eurasian distribution, and theNorth American White Geese (Anser).
Within these three groups, there is substantial variation in beakmorphology which is most pronounced in the Grey Geese of thegenus Anser. The Grey Geese can be regarded as an adaptive radi
ation, because they show ‘‘the evolution of ecological and pheno
typic diversity within a rapidly multiplying lineage” (Schluter,2000). To qualify as an adaptive radiation, four criteria should be
met: (1) common ancestry, (2) rapid speciation, (3) phenotype
environment correlation, and (4) trait utility (Schluter, 2000). The
ﬁrst two criteria (common ancestry and rapid speciation) are supported by our phylogenomic analyses. All members of the Grey
Geese share a common ancestor around 3.4 Mya and diversiﬁed
into several distinct species over a period of about one million
years. A phenotype environment correlation indicates that differ
ences in phenotypes are associated with the use of different
resources or environments. For geese, this correlation has been
conﬁrmed by numerous ﬁeld studies (Carriere et al., 1999; Fox
and Bergersen, 2005; Giroux and Bergeron, 1996). For instance,
the Taiga Bean Goose has a long thin bill adapted for probing soft
substrates of bogs and marshlands to feed on underground plantmaterial, whereas the Pink footed Goose has a short bill suitablefor grazing and seed stripping (Owen, 1976; Reed, 1976). Finally,
the trait utility (i.e., evidence that morphological or physiologicaltraits are indeed useful when used in the associated environment)of beak morphology in geese has been conﬁrmed experimentally(Durant et al., 2003; Heuermann et al., 2011; Zhang et al., 2015).
The Grey Geese can thus be regarded as an adaptive radiation interms of beak morphology.
The diversiﬁcation in beak morphology is less pronounced in
the Black Geese and White Geese, which can be attributed to thefact that some members of the Black Geese (speciﬁcally, theWhite cheeked Geese) and the White Geese coexisted in North
America and consequently affected each other in terms of beakmorphology evolution by means of resource competition and con
sequent character displacement (Dayan and Simberloff, 2005;
Schluter, 2000). So, while the Eurasian niches were rapidly ﬁlled
by members of one ancestral Grey Goose lineage, the North American niches were ﬁlled by members of two distantly related lin
eages: the Black Geese and the White Geese. One can see every
North American species as the ecological equivalent of a Eurasianspecies (or vice versa), for example, Owen (1980) remarked that‘‘the larger races of the Bean Goose [are] the Canada’s [goose] eco
logical counterpart in Eurasia.” However, it is also possible that theBlack Geese and White Geese both exhibited an adaptive radiationin North America comparable with the Eurasian radiation of the
Grey Geese, but that several species went extinct during the Pleis
tocene (Martin and Klein, 1989).
We hypothesized a Pliocene origin for modern goose species
followed by diversiﬁcation at the subspecies level during the Pleis
tocene, based on the available fossil evidence (Brodkorb, 1964) andprevious molecular studies (Avise et al., 1992; Ruokonen et al.,
2000; Scribner et al., 2003; Van Wagner and Baker, 1990). Indeed,
our analysis indicates that the majority of speciation events took
place at the end of the Pliocene. The approximate date of diversiﬁ
cation coincides with the beginning of a period of climatic oscillations between 3.2 and 1.9 million years ago. This period was part ofa fast global cooling trend, following the closure of the Panama
Seaway and the uplifting of the Tibetan Plateau around four million
years ago (Zachos et al., 2001). This resulted in the formation of
permanent Northern Hemisphere ice sheets, the establishment ofa circumpolar tundra belt and the emergence of temperate grasslands, which opened up new ecological niches in which newgroups of animals and plants were able to spread (Kahlke, 2014;
Prins, 1998). The tundra habitat serves as breeding ground forgeese (Owen, 1980), while the temperate grasslands act as wintering grounds where mate choice takes place (Rodway, 2007). Moreover, these tundra and grassland habitats provided ample
opportunity for geese to explore new ecological niches and diversify in beak morphology.
5. Conclusions
Using a phylogenomic approach we were able to resolve thecontentious phylogenetic history of the True Geese. Furthermore,taking advantage of the many contrasting gene histories, we
gained more insight into the effects of complex speciation pro
cesses, such as rapid diversiﬁcation and hybridization, in certain
clades. The relative importance of hybridization in the evolution
ary history of the True Geese remains to be investigated. The wide
spread occurrence of hybridization in birds (Ottenburghs et al.,2015), and speciﬁcally waterfowl (Kraus et al., 2012; Randler,
2008; Ottenburghs et al., 2016), suggests that hybridization can
act as an important component in avian speciation (Rheindt and
Edwards, 2011). By integrating over the full exon set of genes wemade a ﬁrst step to quantitatively describe both species and gene
histories. Our approach will be a fruitful strategy for resolving
many other complex evolutionary histories at the level of genera,
species, and subspecies.
Acknowledgments
We thank three anonymous reviewers for their comments on
the previous version of this manuscript. We are indebted to Gerard
Müskens, Jurre Brenders, Henk en Wim Meinen, Ouwehands Zoo(Gerard Meijer), Avifauna (Jan Harteman & Joost Lammers) andthe Nederlands Instituut voor Ecologie (NIOO KNAW) for their help
during the collection of blood samples, and to Kimberley Laport for
the DNA extraction. We thank Laurent Frantz and Joost de Jong forvaluable discussions and Audrie Siahainenia for the preparation ofFig. 4. This research was funded by Stichting de Eik.
Appendix A. Supplementary material
Supplementary data associated with this article can be found, inthe online version, at http://dx.doi.org/10.1016/j.ympev.2016.05.
021.
References
Ane, C., Larget, B., Baum, D.A., Smith, S.D., Rokas, A., 2007. Bayesian estimation of
concordance among gene trees. Mol. Biol. Evol. 24, 412–426.
Avise, J.C., Alisauskas, R.T., Nelson, W.S., Ankney, C.D., 1992. Matriarchal population
genetic-structure in an avian species with female natal philopatry. Evolution 46,
1084–1096.
Ballard, J.W.O., Whitlock, M.C., 2004. The incomplete natural history of
mitochondria. Mol. Ecol. 13, 729–744.
Baum, D.A., 2007. Concordance trees, concordance factors, and the exploration ofreticulate genealogy. Taxon 56, 417–426.
Berger, A.J., 1972. Hawaiian birdlife.
Bouckaert, R.R., 2010. DensiTree: making sense of sets of phylogenetic trees.
Bioinformatics 26, 1372–1373.
Britton, T., Anderson, C.L., Jacquet, D., Lundqvist, S., Bremer, K., 2007. Estimating
divergence times in large phylogenetic trees. Syst. Biol. 56, 741–752.
Brodkorb, P., 1964. Catalogue of Fossil Birds: Part 2 (Anseriformes throughGalliformes). University of Florida.
Brown, J.W., Rest, J.S., Garcia-Moreno, J., Sorenson, M.D., Mindell, D.P., 2008. Strong
mitochondrial DNA support for a Cretaceous origin of modern avian lineages.
BMC Biol., 6
Bryant, D., 2003. A classiﬁcation of consensus methods for phylogenetics. DIMACS
Series Discr. Math. Theor. Comput. Sci. 61, 163–184.
Carling, M.D., Lovette, I.J., Brumﬁeld, R.T., 2010. Historical divergence and geneﬂow: coalescent analyses of mitochondrial, autosomal and sex-linked loci inpasserina buntings. Evolution 64, 1762–1772.
Carriere, S., Bromley, R.G., Gauthier, G., 1999. Comparative spring habitat and fooduse by two arctic nesting geese. Wilson Bull. 111, 166–180.
Clements, J.F., Schulenberg, T.S., Iliff, M.J., Roberson, D., Fredericks, T.A., Sullivan, B.
L., Wood., C.L., 2015. The eBird/Clements checklist of birds of the world:

312v2015. Downloaded from < http://www.birds.cornell.edu/clementschecklist/
download/>.
DaCosta, J.M., Sorenson, M.D., 2016. DdRAD-seq phylogenetics based on nucleotide,indel, and presence-absence polymorphisms: analyses of two avian genera with
contrasting histories. Mol. Phylogenet. Evol. 94, 122–135.
Davalos, L.M., Cirranello, A.L., Geisler, J.H., Simmons, N.B., 2012. Understanding
phylogenetic incongruence: lessons from phyllostomid bats. Biol. Rev. 87, 991–
1024.
Dayan, T., Simberloff, D., 2005. Ecological and community-wide character
displacement: the next generation. Ecol. Lett. 8, 875–894.
de Queiroz, A., Gatesy, J., 2007. The supermatrix approach to systematics. Trends
Ecol. Evol. 22, 34–41.
Degnan, J.H., Rosenberg, N.A., 2009. Gene tree discordance, phylogenetic inference
and the multispecies coalescent. Trends Ecol. Evol. 24, 332–340.
Delacour, J., 1951. Taxonomic notes on the bean geese, Anser fabalis Lath. Ardea 39,
135–142.
Delacour, J., Mayr, E., 1945. The family Anatidae. Wilson Bull., 3–55
Delsuc, F., Brinkmann, H., Philippe, H., 2005. Phylogenomics and the reconstructionof the tree of life. Nat. Rev. Genet. 6, 361–375.
Donne-Gousse, C., Laudet, V., Hanni, C., 2002. A molecular phylogeny ofanseriformes based on mitochondrial DNA analysis. Mol. Phylogenet. Evol. 23,
339–356.
Drummond, A.J., Suchard, M.A., Xie, D., Rambaut, A., 2012. Bayesian phylogenetics
with BEAUti and the BEAST 1.7. Mol. Biol. Evol. 29, 1969–1973.
Durant, D., Fritz, H., Blais, S., Duncan, P., 2003. The functional response in three
species of herbivorous Anatidae: effects of sward height, body mass and bill
size. J. Anim. Ecol. 72, 220–231.
Ellegren, H., Smeds, L., Burri, R., Olason, P.I., Backström, N., Kawakami, T., Kunstner,
A., Makinen, H., Nadachowska-Brzyska, K., Qvarnström, A., Uebbing, S., Wolf, J.B.
W., 2012. The genomic landscape of species divergence in Ficedula ﬂycatchers.
Advance online publication.
Fox, A.D., Bergersen, L., 2005. Lack of competition between barnacle geese Branta
leucopsis and pink-footed geese Anser brachyrhynchus during the pre-breeding
period in Svalbard. J. Avian Biol. 36, 173–178.
Frantz, L., 2015. Speciation and Domestication in Suiformes: A Genomic Perspective.
Wageningen University, Wageningen.
Frantz, L.A.F., Schraiber, J.G., Madsen, O., Megens, H.J., Bosse, M., Paudel, Y., Semiadi,
G., Meijaard, E., Li, N., Crooijmans, R.P.M.A., Archibald, A.L., Slatkin, M., Schook,
L.B., Larson, G., Groenen, M.A.M., 2013. Genome sequencing reveals ﬁne scale
diversiﬁcation and reticulation history during speciation in Sus. Genome Biol.,14
Fulton, T.L., Letts, B., Shapiro, B., 2012. Multiple losses of ﬂight and recent speciationin steamer ducks. Proc. Roy. Soc. B-Biol. Sci. 279, 2339–2346.
Gil, F., Donsker, D., 2013. IOC World Bird List (v 3.3). Available at <http://
www.worldbirdnames.org>.
Giroux, J.F., Bergeron, R., 1996. Spring diets of sympatric Greater Snow Geese andCanada Geese in southern Quebec. Can. J. Zool. 74, 950–953.
Gonzalez, J., Duttmann, H., Wink, M., 2009. Phylogenetic relationships based on two
mitochondrial genes and hybridization patterns in Anatidae. J. Zool. 279, 310–
318.
Graybeal, A., 1998. Is it better to add taxa or characters to a difﬁcult phylogenetic
problem? Syst. Biol. 47, 9–17.
Heuermann, N., van Langevelde, F., van Wieren, S.E., Prins, H.H., 2011. Increased
searching and handling effort in tall swards lead to a Type IV functional
response in small grazing herbivores. Oecologia 166, 659–669.
Huang, Y.H., Li, Y.R., Burt, D.W., Chen, H.L., Zhang, Y., Qian, W.B., Kim, H., Gan, S.Q.,
Zhao, Y.Q., Li, J.W., Yi, K., Feng, H.P., Zhu, P.Y., Li, B., Liu, Q.Y., Fairley, S., Magor, K.
E., Du, Z.L., Hu, X.X., Goodman, L., Tafer, H., Vignal, A., Lee, T., Kim, K.W., Sheng, Z.
Y., An, Y., Searle, S., Herrero, J., Groenen, M.A.M., Crooijmans, R.P.M.A., Faraut, T.,
Cai, Q.L., Webster, R.G., Aldridge, J.R., Warren, W.C., Bartschat, S., Kehr, S., Marz,
M., Stadler, P.F., Smith, J., Kraus, R.H.S., Zhao, Y.F., Ren, L.M., Fei, J., Morisson, M.,
Kaiser, P., Grifﬁn, D.K., Rao, M., Pitel, F., Wang, J., Li, N., 2013. The duck genome
and transcriptome provide insight into an avian inﬂuenza virus reservoir
species. Nat. Genet. 45, 776–783.
Humphries, E.M., Peters, J.L., Jonsson, J.E., Stone, R., Afton, A.D., Omland, K.E., 2009.
Genetic differentiation between sympatric and allopatric wintering populationsof snow geese. Wilson J. Ornithol. 121, 730–738.
Hung, C.M., Drovetski, S.V., Zink, R.M., 2012. Multilocus coalescence analyses
support a mtDNA-based phylogeographic history for a widespread palearctic
passerine bird, sitta Europaea. Evolution 66, 2850–2864.
Ilves, K.L., Lopez-Fernandez, H., 2014. A targeted next-generation sequencing toolkit
for exon-based cichlid phylogenomics. Mol. Ecol. Resour. 14, 802–811.
Innan, H., Watanabe, H., 2006. The effect of gene ﬂow on the coalescent time in thehuman-chimpanzee ancestral population. Mol. Biol. Evol. 23, 1040–1047.
Jarvis, E.D., Mirarab, S., Aberer, A.J., Li, B., Houde, P., Li, C., Ho, S.Y.W., Faircloth, B.C.,
Nabholz, B., Howard, J.T., Suh, A., Weber, C.C., da Fonseca, R.R., Li, J.W., Zhang, F.,
Li, H., Zhou, L., Narula, N., Liu, L., Ganapathy, G., Boussau, B., Bayzid, M.S.,
Zavidovych, V., Subramanian, S., Gabaldon, T., Capella-Gutierrez, S., HuertaCepas, J., Rekepalli, B., Munch, K., Schierup, M., Lindow, B., Warren, W.C., Ray, D.,
Green, R.E., Bruford, M.W., Zhan, X.J., Dixon, A., Li, S.B., Li, N., Huang, Y.H.,
Derryberry, E.P., Bertelsen, M.F., Sheldon, F.H., Brumﬁeld, R.T., Mello, C.V., Lovell,
P.V., Wirthlin, M., Schneider, M.P.C., Prosdocimi, F., Samaniego, J.A., Velazquez,
A.M.V., Alfaro-Nunez, A., Campos, P.F., Petersen, B., Sicheritz-Ponten, T., Pas, A.,
Bailey, T., Scoﬁeld, P., Bunce, M., Lambert, D.M., Zhou, Q., Perelman, P., Driskell,A.C., Shapiro, B., Xiong, Z.J., Zeng, Y.L., Liu, S.P., Li, Z.Y., Liu, B.H., Wu, K., Xiao, J.,
Yinqi, X., Zheng, Q.M., Zhang, Y., Yang, H.M., Wang, J., Smeds, L., Rheindt, F.E.,

Braun, M., Fjeldsa, J., Orlando, L., Barker, F.K., Jonsson, K.A., Johnson, W., Koepﬂi,
K.P., O’Brien, S., Haussler, D., Ryder, O.A., Rahbek, C., Willerslev, E., Graves, G.R.,
Glenn, T.C., McCormack, J., Burt, D., Ellegren, H., Alstrom, P., Edwards, S.V.,
Stamatakis, A., Mindell, D.P., Cracraft, J., Braun, E.L., Warnow, T., Jun, W., Gilbert,
M.T.P., Zhang, G.J., 2014. Whole-genome analyses resolve early branches in thetree of life of modern birds. Science 346, 1320–1331.
Jetz, W., Thomas, G.H., Joy, J.B., Hartmann, K., Mooers, A.O., 2012. The global
diversity of birds in space and time. Nature 491, 444–448.
Jiang, F.W., Miao, Y.W., Liang, W., Ye, H.Y., Liu, H.L., Liu, B., 2010. The complete
mitochondrial genomes of the whistling duck (Dendrocygna javanica) and blackswan (Cygnus atratus): dating evolutionary divergence in Galloanserae. Mol.
Biol. Rep. 37, 3001–3015.
Jonker, R.M., Kraus, R.H.S., Zhang, Q., Van Hooft, P., Larsson, K., Van der Jeugd, H.P.,
Kurvers, R.H.J.M., Van Wieren, S.E., Loonen, M.J.J.E., Crooijmans, R.P.M.A.,
Ydenberg, R.C., Groenen, M.A.M., Prins, H.H.T., 2013. Genetic consequences ofbreaking migratory traditions in barnacle geese Branta leucopsis. Mol. Ecol. 22,
5835–5847.
Kahlke, R.D., 2014. The origin of Eurasian Mammoth Faunas (MammuthusCoelodonta Faunal Complex). Quat. Sci Rev. 96, 32–49.
Kerr, K.C.R., Cloutier, A., Baker, A.J., 2014. One hundred new universal exonic
markers for birds developed from a genomic pipeline. J. Ornithol. 155, 561–569.
Kingman, J.F.C., 1980. Mathematics of Genetic Diversity. Society for Industrial andApplied Mathematics, Philadelphia, Pa.
Knowles, L.L., 2009. Estimating species trees: methods of phylogenetic analysis
when there is incongruence across genes. Syst. Biol. 58, 463–467.
Knowles, L.L., Kubatko, L.S., 2010. Estimating Species Trees: Practical andTheoretical Aspects. Wiley-Blackwell, Hobroken, NJ.
Kraus, R.H.S., Kerstens, H.H.D., van Hooft, P., Megens, H.J., Elmberg, J., Tsvey, A.,
Sartakov, D., Soloviev, S.A., Crooijmans, R.P.M.A., Groenen, M.A.M., Ydenberg, R.
C., Prins, H.H.T., 2012. Widespread horizontal genomic exchange does not erode
species barriers among sympatric ducks. BMC Evol. Biol., 12
Kubatko, L.S., Carstens, B.C., Knowles, L.L., 2009. STEM: species tree estimation using
maximum likelihood for gene trees under coalescence. Bioinformatics 25, 971–
973.
Lavretsky, P., McCracken, K.G., Peters, J.L., 2014. Phylogenetics of a recent radiationin the mallards and allies (Ayes: Anas): inferences from a genomic transect andthe multispecies coalescent. Mol. Phylogenet. Evol. 70, 402–411.
Lee, S.Y., Scott, G.R., Milsom, W.K., 2008. Have wing morphology or ﬂight kinematics
evolved for extreme high altitude migration in the bar-headed goose? Comp.
Biochem. Phys. C 148, 324–331.
Li, H., Durbin, R., 2009. Fast and accurate short read alignment with BurrowsWheeler transform. Bioinformatics 25, 1754–1760.
Liu, K., Linder, C.R., Warnow, T., 2010. Multiple sequence alignment: a major
challenge to large-scale phylogenetics. PLoS Curr. 2, RRN1198.
Liu, L., Wu, S.Y., Yu, L.L., 2015. Coalescent methods for estimating species trees fromphylogenomic data. J. Syst. Evol. 53, 380–390.
Livezey, B.C., 1996. A phylogenetic analysis of geese and swans (Anseriformes:
Anserinae), including selected fossil species. Syst. Biol. 45, 415–450.
Lu, L., Chen, Y., Wang, Z., Li, X., Chen, W., Tao, Z., Shen, J., Tian, Y., Wang, D., Li, G.,2015. The goose genome sequence leads to insights into the evolution ofwaterfowl and susceptibility to fatty liver. Genome Biol. 16, 89.
Macdonald, G.A., Abbott, A.T., Peterson, F.L., 1983. Volcanoes in the Sea: The
Geology of Hawaii. University of Hawaii Press.
Maddison, W.P., 1997. Gene trees in species trees. Syst. Biol. 46, 523–536.
Martin, P.S., Klein, R.G., 1989. Quaternary Extinctions: A Prehistoric Revolution.
University of Arizona Press.
Mooij, J., Zöckler, C., 1999. Reﬂections on the systematics, distribution and status ofAnser fabalis (Latham, 1787). Casarca 5, 103–120.
Nater, A., Burri, R., Kawakami, T., Smeds, L., Ellegren, H., 2015. Resolving
evolutionary relationships in closely related species with whole-genome
sequencing data. Syst. Biol. 64, 1000–1017.
Ottenburghs, J., Ydenberg, R.C., Van Hooft, P., Van Wieren, S.E., Prins, H.H.T., 2015.
The Avian Hybrids Project: gathering the scientiﬁc literature on avian
hybridization. Ibis 157, 892–894.
Ottenburghs, J., van Hooft, P., van Wieren, S.E., Ydenberg, R.C., Prins, H.H.T., 2016.
Hybridization in geese: a review. Front. Zool. 13, 1–9.
Owen, M., 1976. Factors affecting the distribution of geese in the British Isles.
Wildfowl 27, 143–147.
Owen, M., 1980. Wild Geese of the World: Their Life History and Ecology. Batsford,
London.
Pamilo, P., Nei, M., 1988. Relationships between Gene Trees and Species Trees. Mol.
Biol. Evol. 5, 568–583.
Parham, J.F., Donoghue, P.C.J., Bell, C.J., Calway, T.D., Head, J.J., Holroyd, P.A., Inoue, J.
G., Irmis, R.B., Joyce, W.G., Ksepka, D.T., Patane, J.S.L., Smith, N.D., Tarver, J.E., vanTuinen, M., Yang, Z.H., Angielczyk, K.D., Greenwood, J.M., Hipsley, C.A., Jacobs, L.,
Makovicky, P.J., Muller, J., Smith, K.T., Theodor, J.M., Warnock, R.C.M., Benton, M.
J., 2012. Best practices for justifying fossil calibrations. Syst. Biol. 61, 346–359.
Paxinos, E.E., James, H.F., Olson, S.L., Sorenson, M.D., Jackson, J., Fleischer, R.C., 2002.
MtDNA from fossils reveals a radiation of Hawaiian geese recently derived fromthe Canada goose (Branta canadensis). Proc. Natl. Acad. Sci. U. S. A. 99, 1399–
1404.
Payseur, B.A., 2010. Using differential introgression in hybrid zones to identify
genomic regions involved in speciation. Mol. Ecol. Resour. 10, 806–820.
Pereira, S.L., Baker, A.J., 2006. A mitogenomic timescale for birds detects variablephylogenetic rates of molecular evolution and refutes the standard molecular
clock. Mol. Biol. Evol. 23, 1731–1740.

313
Prins, H.H.T., 1998. The origins and development of grassland communities innorthwestern Europe. In: Wallis de Vries, M.F., Bakker, J.P., van Wieren, S.E.
(Eds.), Grazing and Conservation Management. Kluwer Academic Publishers,
Boston, pp. 55–105.
Prum, R.O., Berv, J.S., Dornburg, A., Field, D.J., Townsend, J.P., Lemmon, E.M.,
Lemmon, A.R., 2015. A Comprehensive Phylogeny of Birds (Aves) Using Targeted
Next-Generation DNA Sequencing. Advance online publication.
Quinn, T., 1992. The genetic legacy of Mother Goose–phylogeographic patterns oflesser snow goose Chen caerulescens caerulescens maternal lineages. Mol. Ecol.
1, 105–117.
Quinn, T.W., Shields, G.F., Wilson, A.C., 1991. Afﬁnities of the Hawaiian Goose Basedon 2 Types of Mitochondrial-DNA Data. Auk 108, 585–593.
Randler, C., 2008. Hybrid wildfowl in Central Europe – an overview. Waterbirds 31,
143–146.
Rannala, B., Yang, Z.H., 2003. Bayes estimation of species divergence times andancestral population sizes using DNA sequences from multiple loci. Genetics
164, 1645–1656.
Ranwez, V., Delsuc, F., Ranwez, S., Belkhir, K., Tilak, M.K., Douzery, E.J.P., 2007.
OrthoMaM: a database of orthologous genomic markers for placental mammal
phylogenetics. BMC Evol. Biol., 7
Reed, A., 1976. Geese, nutrition and farmland. Wildfowl 27, 153–155.
Rheindt, F.E., Edwards, S.V., 2011. Genetic introgression: an integral but neglectedcomponent of speciation in birds. Auk 128, 620–632.
Rodway, M.S., 2007. Timing of pairing in waterfowl I: Reviewing the data andextending the theory. Waterbirds 30, 488–505.
Rokas, A., Williams, B.L., King, N., Carroll, S.B., 2003. Genome-scale approaches toresolving incongruence in molecular phylogenies. Nature 425, 798–804.
Ruokonen, M., Aarvak, T., Madsen, J., 2005. Colonization history of the high-arctic
pink-footed goose Anser brachyrhynchus. Mol. Ecol. 14, 171–178.
Ruokonen, M., Kvist, L., Lumme, J., 2000. Close relatedness between mitochondrial
DNA from seven Anser goose species. J. Evol. Biol. 13, 532–540.
Ruokonen, M., Litvin, K., Aarvak, T., 2008. Taxonomy of the bean goose-pink-footed
goose. Mol. Phylogenet. Evol. 48, 554–562.
Sangster, G., Oreel, G.J., 1996. Progress in taxonomy of taiga and tundra bean geese.
Dutch Birding 18, 310–316.
Schluter, D., 2000. The Ecology of Adaptive Radiation. Oxford University Press.
Scribner, K.T., Talbot, S.L., Pearce, J.M., Pierson, B.J., Bollinger, K.S., Derksen, D.V.,
2003. Phylogeography of Canada Geese (Branta canadensis) in western North
America. Auk 120, 889–907.

Shields, G.F., 1990. Analysis of Mitochondrial-DNA of Paciﬁc Black Brant (BrantaBernicla-Nigricans). Auk 107, 620–623.
Shields, G.F., Wilson, A.C., 1987. Calibration of Mitochondrial-DNA Evolution in
Geese. J. Mol. Evol. 24, 212–217.
Smith, N.A., 2015. Sixteen vetted fossil calibrations for divergence dating of
Charadriiformes (Aves, Neognathae). Palaeontol Electr., 18
Stamatakis, A., 2006. RAxML-VI-HPC: maximum likelihood-based phylogenetic
analyses with thousands of taxa and mixed models. Bioinformatics 22, 2688–
2690.
Toews, D.P.L., Brelsford, A., 2012. The biogeography of mitochondrial and nuclear
discordance in animals. Mol. Ecol. 21, 3907–3930.
van Tuinen, M., Hedges, S.B., 2001. Calibration of avian molecular clocks. Mol. Biol.
Evol. 18, 206–213.
Van Wagner, C.E., Baker, A.J., 1990. Association between Mitochondrial-DNA andMorphological Evolution in Canada Geese. J. Mol. Evol. 31, 373–382.
Volkovsky, D., Kulikova, I., Gerasimov, Y.N., Zhuravlev, Y.N., 2013. Genetic diversityof Anser albifrons Scopoli, 1769 and Anser fabalis Latham, 1787 in the Russian
Far East. Russ. J. Genet+ 49, 428–440.
Wall, J.D., 2003. Estimating ancestral population sizes and divergence times.
Genetics 163, 395–404.
Weir, J.T., Schluter, D., 2008. Calibrating the avian molecular clock. Mol. Ecol. 17,
2321–2328.
Wetmore, A., 1956. A Check-List of the Fossil and Prehistoric Birds of North America
and the West Indies. Smithsonian Institution.
Wu, Y.F., 2012. Coalescent-based species tree inference from gene tree topologies
under incomplete lineage sorting by maximum likelihood. Evolution 66, 763–
775.
Yang, Z., Goldman, N., Friday, A., 1994. Comparison of models for nucleotide
substitution used in maximum-likelihood phylogenetic estimation. Mol. Biol.
Evol. 11, 316–324.
Yang, Z.H., 2007. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol.
Evol. 24, 1586–1591.
Zachos, J., Pagani, M., Sloan, L., Thomas, E., Billups, K., 2001. Trends, rhythms, andaberrations in global climate 65 Ma to present. Science 292, 686–693.
Zhang, Y., Jia, Q., Prins, H.H., Cao, L., de Boer, W.F., 2015. Individual-Area
Relationship Best Explains Goose Species Density in Wetlands.