BACKGROUND
flamingos are a unique order  of birds with a highly specialized ecology but their evolutionary history remains poorly understood and until recently has only been informed by the fossil record.
flamingos are traditionally perceived as among the oldest lineages of living birds with reports of flamingo-like birds appearing in the fossil record as early as the late cretaceous .
however, the earliest birds reliably placed as stem phoenicopteriforms  first appear in the early oligocene of europe  <cit>  and the earliest members of the crown family  appear during the oligo-miocene of the old world  <cit> , suggesting an age for the family on par with most other major familial divergences within aves  <cit> .
notably, harrisonavis croizeti, an apparently morphologically modern flamingo from the oligo-miocene of france, suggests the modern flamingo divergence occurred in the old world as early as the late paleogene  <cit> .
recent fossil and molecular work have shed new insight into phoenicopterid origins and cast doubt on the classical notion of flamingos as a particularly ancient lineage among the storks, herons and ibises .
most notably, molecular  <cit>  and morphological  <cit>  studies have supported a sister relationship between flamingos and grebes  as the clade mirandornithes  <cit> .
however, despite these advances, the exact age and phylogeny of modern flamingos remains to be robustly tested and several questions about flamingo evolution are in need of further investigation:  how do the six extant species of flamingos relate to each other,  how long ago did these divergences occur, and  where did crown phoenicopteridae originate?
sibley and ahlquist  <cit>  represents the only genetic study that addresses these questions.
this study, based on differences in hybridization strength between genomic dna of different species, documented a shallow age among five species and  two subspecies  that likely form a single genus comprising two distinct sub-clades.
this division is congruent with the organization by jenkin  <cit>  of flamingos into two groups reflecting mandibular morphology and feeding strategy: the lesser, andean and james’s flamingos  have bulbous bills in cross-section suited to filtration of smaller food items ; the remaining species  have more compact bills in cross-section suited to filtration of larger food items .
further morphological variation within the two sub-clades is not completely known  <cit> .
intraspecific genomic variation in each flamingo species is also incompletely known but is an important consideration when investigating morphological and geographic origins of species with wide ranges.
recent population genetics studies have identified a lack of genetic structure within the greater  <cit> ; remi wattier, pers.
comm.
and lesser  <cit>  flamingos, both species with extensive ranges in the old world.
it is unlikely that the restricted breeding distribution  <cit>  on the andean altiplano promotes genetic structure in the andean and james’s flamingos.
except for captive populations  <cit> , the population genetics of the american and chilean flamingos remains uninvestigated.
the present study has three aims: 1) to assess intraspecific mitochondrial variation within the widespread new world chilean flamingo through the use of three mtdna markers and museum skins collected throughout the range; 2) to test the phylogenetic relationships among modern flamingos through the use of a multi-locus  and varied analytical  approach; and 3) to estimate the divergence times of living flamingo species from the resulting species tree calibrated with two fossil constraints.
how these results inform the taxonomic organization and biogeography of flamingos is discussed.
RESULTS
population genetics of the chilean flamingo
sequence data for  <dig> mtdna loci were obtained and analyzed from 17– <dig> chilean flamingos representing the entire breeding range of the species except ecuador .
only two out of  <dig> individuals showed variation in the coi gene , exposing the same substitution recovered previously by a dna barcoding study of neotropical birds  <cit> .
that study indicates the presence of only one substitution  in the entire barcode portion of the coi.
from the current study, the common haplotype was distributed throughout the range, while the minor variant was found in an individual from peru  and in three individuals from the argentinean andes .
both haplotypes coexist in these localities.
two additional markers that are known to be variable in related species  failed to reveal the existence of additional mitochondrial variation.
phylogenetics
analysis of single loci generally failed to result in complete phylogenetic resolution, instead supporting one or few of the key nodes in the final species tree .
none of the single locus-based gene trees  indicate significant conflict, rather lack of resolution, with the final multi locus-based species tree.
thus, the joint analyses of loci that vary in mutation rate contributed to enhanced overall resolution the species topology is congruent with previous non-dna sequence based studies  <cit>  with flamingos breaking into two distinct clades: a shallow-keeled clade comprising ) and a deep-keeled clade comprising ).
additional analysis indicates that the species topology is particularly robust.
the bayes factor analysis supported partitioning of the total dataset by locus , though both partitioned and non-partitioned bi analyses recovered the same topology with complete support ; the ml analysis had >85% bootstrap support for all nodes.
a *beast analysis including multiple individuals for each flamingo species again recovered the same topology with complete support .
all jackknife analyses excluding a single locus at a time from the total  dataset recovered the complete topology with high  posterior probability, suggesting that no single locus was driving the results.
the ndna dataset recovered the same topology, while the mtdna dataset differed only by recovering a ) clade, albeit without significant support.
even though intra-specific variation was discovered in a set of five ndna loci including two to three individuals per flamingo species, this variation was found to be minimal , an inference confirmed by independent cross-checking between the two genetics labs which revealed a match for all  <dig> flamingos.
furthermore, five of the  <dig> loci showed interspecific variation between the american and greater flamingos , with the barcode locus  k2p distance indicating 1% variation.
divergence time estimation
the dates recovered by the mtdna dataset are consistently ≥ <dig> times the ages recovered by the ndna dataset, while those recovered by the total dataset fall in between .
crown-clade phoenicopteridae likely arose at or shortly after the mio-pliocene boundary .
the deep- and shallow-keeled clades diverged in either the pliocene or earliest pleistocene  followed by the american-greater split  and the andean-james’s split .
age ranges represent the total  and the ndna only  datasets.
see table  <dig> for divergence dating results of the all datasets .
ages of divergences within flamingos are in bold followed by the 95% confidence intervals in parentheses.
extant flamingos likely originated at or shortly after the miocene-pliocene boundary, making them one of the youngest lineages among living birds.
the ages recovered by the nuclear data are consistently much younger than the ages recovered by the mitochondrial data, likely the result of mutational saturation in the mitochondrial sequences.
discussion
population genetics of the chilean flamingo
a single metapopulation was identified across the range of the chilean flamingo, harboring minimal genetic variation.
this finding parallels similar recent findings in other widespread flamingo species  <cit>  and is consistent with the appearance of opportunistic breeding  <cit>  and a nomadic lifestyle  <cit>  in all flamingos.
given these results, minimal genetic structure is expected for the phylogeographically unstudied andean, james’s and american flamingos.
thus, it is unlikely that unaccounted-for population structure impacted the results of our phylogenetic and divergence dating analyses.
phylogenetics
our study provides the first sequence-based molecular support for the relationships among flamingos suggested previously  <cit> .
by contrasting the lack of full resolution for any of the individual loci with the complete resolution found for the ndna and total datasets, the need for multiple independent loci to resolve temporally young divergences is underscored.
conflicting phylogenetic signal exists between the mtdna loci as coi was the only single locus to fully resolve the flamingo relationships while cyt b recovered the andean flamingo at the base of the deep-keeled clade with high support.
the relationships recovered by the total dataset agreed with those recovered by coi, as well as by the ndna dataset and those recovered previously  <cit> , indicating that the ndna data sufficiently overpowered the conflicting cyt b signal.
the two sub-clades recovered in our analysis are congruent with the patterns detected by jenkin  <cit>  in mandibular morphology.
the shallow-keeled clade comprises those species with bills suited to capture larger prey items and includes the greater , american  and chilean  flamingos.
within this clade, the american and greater flamingos are often considered sub-species .
the 1% divergence inferred from the dna barcode locus falls in the range of expected nearest neighbor distances among sister species, such as seen commonly in waterfowl  <cit> .
although more sampling is needed to confirm species status, our finding of fixed mutational differences in several nuclear loci and a 1% divergence in mtdna is consistent with treating the american and greater flamingos as separate species.
the deep-keeled clade comprises those species with bills suited to capturing smaller prey items and includes the lesser , andean  and james’s  flamingos.
these species have often been set apart from the classical genus phoenicopterus with the lesser placed within phoeniconaias based on unique characteristics of the bill  <cit>  and the andean and james’s placed together in phoenicoparrus based on the shared absence of the hallux   <cit> .
however, the morphological definition of phoeniconaias does not seem to exclude andinus and jamesi, and the hallux of minor is reduced  <cit> , approximating the morphological definition of phoenicoparrus.
thus, we suggest that the division of the deep-keeled clade into two genera is arbitrary and uninformative.
instead, we suggest the placement of minor within phoenicoparrus along with andinus and jamesi based on shared mandibular morphology, ecology and phylogenetic relatedness.
we suggest the designation of phoeniconaias as a junior synonym of phoenicoparrus based on priority and the redefinition of phoenicoparrus as all species more closely related to phoenicoparrus andinus than to phoenicopterus roseus.
divergence time estimation
even though the various subsets of data recovered incongruent divergence times , all agreed on a recent divergence for crown phoenicopteridae, indicating that living flamingos are among the youngest lineages within neoaves.
this young age is counter to the general perception of the phoenicopteridae as among the oldest lineages of living birds based on fossil data  but is congruent with the dna-dna hybridization data of sibley and ahlquist  <cit> , which identified the close genetic relatedness of the flamingos and predicted an age of 5– <dig> ma for the basal divergence among the extant species.
the ndna dataset recovered ages significantly  younger than those of the mtdna dataset.
this pattern is likely due in part to artificial signal caused by saturation of the more quickly mutating mtdna genes  <cit> .
though saturation is unlikely to be an issue with divergences as shallow as those of crown phoenicopteridae, the use of temporally more distant calibrations external to the clade of interest  may have introduced saturation artifact.
it is likely that homoplasy within the relatively longer and more quickly-evolving mtdna sequences, with more mutational sites, conflicted with the signal of the shorter and more slowly evolving ndna sequences, resulting in artificially old age estimates.
thus, the true ages of the component flamingo divergences are probably somewhere between those predicted by the total dataset and by the ndna dataset alone.
biogeography
the two old world species  are each recovered in clades with otherwise exclusively new world distributions.
the recovery of the greater flamingo as derived within phoenicopterus suggests a new world origin for the shallow-keeled flamingos.
the origin of phoenicoparrus is not so straightforward, as the basal divergence within the clade spans the atlantic.
plio-pleistocene fossil flamingos are known only from the new world  <cit>  and australia  <cit>  and the phylogenetic affinities of these species are uncertain.
no flamingos occupy australia today and it is unknown if the plio-pleistocene fossils represent an extension of the current distribution or if they are descendants of the species which occupied australia in the oligo-miocene  <cit> .
the new world fossil record, however, is congruent with the modern flamingo distribution and likely suggests an origin for phoenicoparrus, and crown phoenicopteridae, in the western hemisphere, followed by trans-atlantic dispersal.
better phylogenetic resolution for the australian forms will further inform phoenicopterid origins.
implications for the fossil record
the phoenicopteriform fossil record is rich but poorly understood and revision will greatly advise divergence time estimation and paleo-biogeographical reconstruction for crown flamingos.
the stem phoenicopteriform record is characterized by fragmentary material with only tentative assignment to the flamingos.
while the earliest palaelodids  are the earliest to be reliably placed in phoenicopteriformes, juncitarsus <cit>  , elornis and agnopterus <cit>  from the middle to late eocene of north america and europe may be earlier members and perhaps more appropriate fossil calibrations for mirandornithes.
better understanding of the plio-pleistocene record may also provide internal calibrations for crown flamingos.
finally, while we present the first information on the timing of diversification and specialization within crown flamingos, better understanding of interspecific variation  is necessary to interpret the ecological significance of these results.
likewise, revision of the stem phoenicopteriform fossil record is necessary to place these ages in a more meaningful context.
most notably, the primitive flamingo harrisonavis croizeti from the oligo-miocene of europe displayed a mandibular morphology approximating that of modern flamingos  <cit> .
reassessment of this material, and new descriptions of contemporaneous taxa, will inform the state of flamingo specialization at this point and better constrain the rate of specialization within the phoenicopteridae.
CONCLUSIONS
chilean flamingos show no evidence of population structure across their entire range, a condition similar to other widely-distributed flamingo species.
the six living species of flamingos fall into two clades reflecting differences in bill morphology and feeding ecology.
we propose slight reorganization of flamingo taxonomy to indicate these differences: the greater, american and chilean flamingos comprise genus phoenicopterus while the lesser is moved to genus phoenicoparrus  along with the james’s and andean flamingos.
the living species diverged from each other ~ <dig>  million years ago, followed by phoenicoparrus  and phoenicopterus .
methods and materials
population genetics methods
field sampling of all flamingo populations remains impractical for logistic reasons, including geographic inaccessibility and conservation protection.
sampling from captive birds alone is not guaranteed to capture species-wide variation due to a persistent lack of information regarding the geographic origins of the ancestors of captive populations.
to assess phylogeographic structure, we thus set out to sample museum specimens of chilean flamingos with known locality.
available dna barcodes from chilean flamingos suggest the presence of at least two haplotypes in argentina  <cit>  but no data is available from other portions of their range.
thus, we examined portions of three mitochondrial genes through amplification from historic toe pads from chilean flamingos  with known localities .
the samples used spanned nearly the entire range known of chilean flamingos  and were obtained through several museum loans.
the three loci were chosen based on known intraspecific variation in flamingos  or known high rates of mutation .
primers were manually generated using sequences available in genbank .
in order to prevent contamination, ancient dna procedures for extraction and amplification were followed  <cit>  .
taxonomic sampling
tissue samples were obtained for each flamingo species .
podilymbus podiceps and aechmophorus clarkii were included to represent the basal split within grebes  <cit> .
the sister taxon to  remains controversial but the tropicbirds  are consistently found as close relatives  and were used as an outgroup.
penguins  and tubenoses  were included as representatives of the waterbird clade .
sample collection was reviewed and approved by the centre de recherche sur la biologie des populations d’oiseaux  under the personal permit  of alan johnson and arnaud béchet and by usda permit  of mvt.
genomic sampling
genomic dna was extracted from tissue samples using the dneasy tissue kit  and was amplified for  <dig> nuclear and two mitochondrial loci from one individual per species, which comprised the primary dataset .
two additional nuclear loci were also amplified to test intraspecific variation among flamingos only.
see additional file  <dig> for primer sequences.
the nuclear loci consisted primarily of intronic regions with primers designed using an exon-primed intron spanning approach and sequences from the chicken and zebra finch genome .
polymerase chain-reaction  amplifications were carried out in  <dig> µl reactions comprising:  <dig> µl dh2o,  <dig> µl 10× detergent- and mg+-free reaction buffer,  <dig> µl mgcl <dig> ,  <dig> µl dntp ,  <dig>  µl each primer,  <dig>  µl hot multitaq  polymerase  and  <dig>  µl dna template .
thermocycling comprised an initial denaturation at 95°c for  <dig> min followed by 35– <dig> cycles of annealing at 95°c for  <dig> sec, extension at 50°c for  <dig> sec and denaturation at 95°c for  <dig> sec, followed by a final annealing/extension at 72°c for  <dig> min.
excess dntps and primers were removed from pcr product by addition of  <dig>  µl exonuclease  <dig> ,  <dig>  µl fastap  shrimp alkaline phosphatase  and  <dig>  µl dh <dig> and held at 37°c for  <dig> min and then 80°c for  <dig> min.
product was sequenced offsite at macrogen usa .
only the two mtdna loci were treated as coding; the mtdna were further partitioned by codon in the divergence time estimation analyses.
substitution models for each locus were estimated using mega  <dig>  and the best available model was used in mega  <dig>  and mrbayes  <dig> .
<dig> for the phylogenetic analyses and in beast  <dig> .
<dig> for the divergence time estimation analyses.
adamts <dig> and hmgb <dig> were only used in phylogenetic analysis testing flamingo intraspecific variation.
abbreviations: rheb <dig>  rheb isoform  <dig>  intron 3; timm17a, mitochondrial import inner membrane translocase subunit tim17a, intron 3; tcf <dig>  transcription factor  <dig>  intron 12; rps <dig>  ribosomal protein s <dig>  intron 5; slc29a <dig>  solute carrier family  <dig>  member  <dig>  intron 8; nfkbiz, nuclear factor of light polypeptide gene enhancer in b-cells inhibitor, zeta, intron 6; g3pdh, glyceraldehyde-3-phosphate dehydrogenase gene, intron 11; myoglobin gene, exon  <dig>   <dig> and intron 2; zenk, zinc finger protein, exon 2; zenk 3’utr, zinc finger protein, 3’ untranslated region; coi, cytochrome oxidase subunit i; cyt b, cytochrome b; adamts <dig>  adam metallopeptidase with thrombospondin, type  <dig> motif, intron 5; hmgb <dig>  high mobility protein group, box  <dig>  intron  <dig>  drawings by m. mccracken.
sequence data for one to two additional flamingo individuals per species for five loci  were collected in the dijon, france lab.
for these samples, genomic dna was extracted from blood samples for the greater, american and chilean flamingos and feather samples for the remaining three species using a standard phenol-chloroform method.
pcr amplifications and product clean-up followed the same protocols as the wilmington, us lab.
product was sequenced offsite at macrogen korea .
sequence treatment and phylogenetic analysis
sequences were aligned within sequencher  <dig>  and alignments were confirmed within mega  <dig>   <cit>  using clustalw with default settings.
pairwise genetic distances were calculated in mega  <dig>  using the default maximum composite likelihood substitution model except for coi, which was calculated using kimura 2-parameter.
maximum likelihood phylogenetic analyses were carried out within mega  <dig>  and the loci were concatenated.
unless the dataset comprised only coding genes , each dataset was treated as noncoding; optimal substitution models were found within the hierarchical bi and ml framework in mega  <dig>   and the analysis used a nearest-neighbor-interchange heuristic search.
bootstrap support was estimated from  <dig> replicates using the same settings.
primary phylogenetic analyses using bayesian inference were carried out in mrbayes  <dig> .
<dig>  <cit> .
datasets were fully partitioned based on the optimal substitution models estimated within mega  <dig> .
each search comprised two concurrent runs of four chains each for  <dig> , <dig> generations sampled every  <dig>  generations with the first  <dig> , <dig> generations discarded as burn-in.
the results of each analysis were tested for convergence of phylogenetic signal in tracer  <dig>   <cit> , where estimated sample size  values > <dig> were treated as reliable signal.
the primary analysis employed all  <dig> ndna and two mtdna loci for one individual per species.
all secondary analyses employed the same parameters as the primary analysis except where noted.
to test for the presence of intraspecific variation, a series of single-locus analyses were carried out using two to three sequences per flamingo species for five loci .
to recover individual gene trees, each locus was analyzed independently.
to test the relative weight of each locus to the combined dataset, a series of analyses were run excluding a single gene at a time.
to test for conflicting signal between the ndna and mtdna datasets, each was analyzed independently.
to test for the effects of over-partitioning of the bayesian inference analysis, a bayesian analysis was carried out treating the total dataset as non-partitioned .
a bayes factor analysis was carried out in tracer  <dig>  to compare the strength in phylogenetic signal of the fully partitioned and non-partitioned bi treatments.
*beast  <dig> .
<dig>  <cit>  was used to test for discordance between individual gene trees and the species tree.
the dataset was fully phased by individual and gene.
loci were allowed unlinked substitution models, clock models and partition trees.
the optimal available substitution model was used  and the  <dig> coding loci were partitioned by codon position.
each locus was allowed a lognormal relaxed clock with uncorrelated rates.
mean clock rate priors  were allowed uninformative uniform distributions .
a yule speciation process was assumed.
no fossil node calibrations were used.
the results were tested for convergence of phylogenetic signal in tracer  <dig> .
divergence timing analysis and fossil calibrations
divergence time estimation analyses were carried out using beast  <dig> .
<dig>  <cit> .
loci were allowed unlinked substitution models  and the  <dig> coding loci were partitioned by codon position.
a lognormal relaxed clock with uncorrelated rates was assumed and a yule speciation process was employed.
the clock mean  was allowed an uninformative uniform prior distribution .
taxa were divided into sets representing grebes and mirandornithes; monophyly of these sets was not enforced.
two nodes were fossil calibrated following the guidelines suggested by parham et al.
<cit> .
the divergence of extant grebes was calibrated at  <dig>  mya based on thiornis sociata, a grebe from the late miocene of spain  <cit> .
the divergence of flamingos and grebes was calibrated at  <dig>  mya based on adelalopus hoogbutseliensis, a stem phoenicopteriform from the earliest oligocene of belgium  <cit> .
calibrations were allowed lognormal distributions with offset = calibration age and standard deviation =  <dig> .
see table  <dig> for the data included in each divergence analysis.
all analyses were run twice for  <dig> million generations sampled every  <dig> with burn-in of  <dig> million.
convergence of phylogenetic signal was tested for in tracer  <dig> .
availability of supporting data
all supporting data are included as additional files.
newly obtained dna sequence reads have been uploaded to genbank .
genbank accession numbers are provided in additional files  <dig> and  <dig>
competing interests
the authors declare that they have no competing interests.
authors’ contributions
crt carried out dna extraction, amplification, alignment, performed all phylogenetic and divergence time analysis, and drafted the manuscript.
lmo and mafg carried out dna extraction, amplification, and alignment, and helped draft the manuscript.
bf carried out ancient dna extraction, amplification and alignment.
mvt conceived of the study, participated in its design and coordination, assisted with alignment and analysis, and helped to draft the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
specimen and locality information for the chilean flamingo individuals included in the population genetics analysis.
click here for file
additional file 2
pairwise genetic distances between multiple individuals of each flamingo species for five ndna loci.
click here for file
additional file 3
pairwise genetic distances between
phoenicopterus ruber
and
p. roseus
for  <dig> ndna and two mtdna loci.
click here for file
additional file 4
primer sequences for each locus used in this study.
click here for file
additional file 5
ancient dna extraction and amplification methods.
click here for file
additional file 6
sample information for individuals used in the primary phylogenetic analyses.
click here for file
additional file 7
genbank accession numbers for sequence data used in the primary phylogenetic analyses.
click here for file
additional file 8
sample information for the additional flamingo individuals used in the pairwise genetic distance analyses and genbank accession numbers for associated sequence data.
click here for file
acknowledgements
we would like to thank the many people and institutes for sample collection and/or approving sample loans: sharon birks , brooks childress and rebecca lee , james dean , jack demboski , becky desjardins and john gerwin , steve emslie , maureen flannery and jack dumbacher , blair hedges , janet hinshaw , peter houde , diana humple , irby lovette , steve cardiff, donna dittman and van remsen , mark robbins , paul sweet and joel cracraft , krystof ziskowski , sylvan heights waterfowl centre , parc des oiseaux , parc zoologique d’amnéville , salins group and tour du valat .
we further like to thank frank cézilly, arnaud béchet and rémi wattier for access to samples and facilities at the université de bourgogne, dijon, france, and dan ksepka for discussion of the flamingo and grebe fossil record.
for financial support, we thank the college of arts and sciences and department of biology and marine biology at uncw , the conseil régional de bourgogne  and the total foundation .
we thank m. mccracken for providing flamingo drawings.

 flamingos are a unique order of birds with a highly specialized ecology but their evolutionary history remains poorly understood and until recently has only been informed by the fossil record. flamingos are traditionally perceived as among the oldest lineages of living birds with reports of flamingo-like birds appearing in the fossil record as early as the late cretaceous. however, the earliest birds reliably placed as stem phoenicopteriforms first appear in the early oligocene of europe and the earliest members of the crown family appear during the oligo-miocene of the old world, suggesting an age for the family on par with most other major familial divergences within aves. notably, harrisonavis croizeti, an apparently morphologically modern flamingo from the oligo-miocene of france, suggests the modern flamingo divergence occurred in the old world as early as the late paleogene. recent fossil and molecular work have shed new insight into phoenicopterid origins and cast doubt on the classical notion of flamingos as a particularly ancient lineage among the storks, herons and ibises. most notably, molecular and morphological studies have supported a sister relationship between flamingos and grebes as the clade mirandornithes. however, despite these advances, the exact age and phylogeny of modern flamingos remains to be robustly tested and several questions about flamingo evolution are in need of further investigation: how do the six extant species of flamingos relate to each other, how long ago did these divergences occur, and where did crown phoenicopteridae originate? sibley and ahlquist represents the only genetic study that addresses these questions. this study, based on differences in hybridization strength between genomic dna of different species, documented a shallow age among five species and two subspecies that likely form a single genus comprising two distinct sub-clades. this division is congruent with the organization by jenkin of flamingos into two groups reflecting mandibular morphology and feeding strategy: the lesser, andean and james’s flamingos have bulbous bills in cross-section suited to filtration of smaller food items; the remaining species have more compact bills in cross-section suited to filtration of larger food items. further morphological variation within the two sub-clades is not completely known. intraspecific genomic variation in each flamingo species is also incompletely known but is an important consideration when investigating morphological and geographic origins of species with wide ranges. recent population genetics studies have identified a lack of genetic structure within the greater; remi wattier, pers. comm. and lesser flamingos, both species with extensive ranges in the old world. it is unlikely that the restricted breeding distribution on the andean altiplano promotes genetic structure in the andean and james’s flamingos. except for captive populations, the population genetics of the american and chilean flamingos remains uninvestigated. the present study has three aims: 1) to assess intraspecific mitochondrial variation within the widespread new world chilean flamingo through the use of three mtdna markers and museum skins collected throughout the range; 2) to test the phylogenetic relationships among modern flamingos through the use of a multi-locus and varied analytical approach; and 3) to estimate the divergence times of living flamingo species from the resulting species tree calibrated with two fossil constraints. how these results inform the taxonomic organization and biogeography of flamingos is discussed.  population genetics of the chilean flamingo sequence data for mtdna loci were obtained and analyzed from 17– chilean flamingos representing the entire breeding range of the species except ecuador. only two out of individuals showed variation in the coi gene, exposing the same substitution recovered previously by a dna barcoding study of neotropical birds. that study indicates the presence of only one substitution in the entire barcode portion of the coi. from the current study, the common haplotype was distributed throughout the range, while the minor variant was found in an individual from peru and in three individuals from the argentinean andes. both haplotypes coexist in these localities. two additional markers that are known to be variable in related species failed to reveal the existence of additional mitochondrial variation. phylogenetics analysis of single loci generally failed to result in complete phylogenetic resolution, instead supporting one or few of the key nodes in the final species tree. none of the single locus-based gene trees indicate significant conflict, rather lack of resolution, with the final multi locus-based species tree. thus, the joint analyses of loci that vary in mutation rate contributed to enhanced overall resolution the species topology is congruent with previous non-dna sequence based studies with flamingos breaking into two distinct clades: a shallow-keeled clade comprising ) and a deep-keeled clade comprising ). additional analysis indicates that the species topology is particularly robust. the bayes factor analysis supported partitioning of the total dataset by locus, though both partitioned and non-partitioned bi analyses recovered the same topology with complete support; the ml analysis had >85% bootstrap support for all nodes. a *beast analysis including multiple individuals for each flamingo species again recovered the same topology with complete support. all jackknife analyses excluding a single locus at a time from the total dataset recovered the complete topology with high posterior probability, suggesting that no single locus was driving the results. the ndna dataset recovered the same topology, while the mtdna dataset differed only by recovering a ) clade, albeit without significant support. even though intra-specific variation was discovered in a set of five ndna loci including two to three individuals per flamingo species, this variation was found to be minimal, an inference confirmed by independent cross-checking between the two genetics labs which revealed a match for all flamingos. furthermore, five of the loci showed interspecific variation between the american and greater flamingos, with the barcode locus k2p distance indicating 1% variation. divergence time estimation the dates recovered by the mtdna dataset are consistently ≥ times the ages recovered by the ndna dataset, while those recovered by the total dataset fall in between. crown-clade phoenicopteridae likely arose at or shortly after the mio-pliocene boundary. the deep- and shallow-keeled clades diverged in either the pliocene or earliest pleistocene followed by the american-greater split and the andean-james’s split. age ranges represent the total and the ndna only datasets. see table for divergence dating results of the all datasets. ages of divergences within flamingos are in bold followed by the 95% confidence intervals in parentheses. extant flamingos likely originated at or shortly after the miocene-pliocene boundary, making them one of the youngest lineages among living birds. the ages recovered by the nuclear data are consistently much younger than the ages recovered by the mitochondrial data, likely the result of mutational saturation in the mitochondrial sequences. discussion population genetics of the chilean flamingo a single metapopulation was identified across the range of the chilean flamingo, harboring minimal genetic variation. this finding parallels similar recent findings in other widespread flamingo species and is consistent with the appearance of opportunistic breeding and a nomadic lifestyle in all flamingos. given these results, minimal genetic structure is expected for the phylogeographically unstudied andean, james’s and american flamingos. thus, it is unlikely that unaccounted-for population structure impacted the results of our phylogenetic and divergence dating analyses. phylogenetics our study provides the first sequence-based molecular support for the relationships among flamingos suggested previously. by contrasting the lack of full resolution for any of the individual loci with the complete resolution found for the ndna and total datasets, the need for multiple independent loci to resolve temporally young divergences is underscored. conflicting phylogenetic signal exists between the mtdna loci as coi was the only single locus to fully resolve the flamingo relationships while cyt b recovered the andean flamingo at the base of the deep-keeled clade with high support. the relationships recovered by the total dataset agreed with those recovered by coi, as well as by the ndna dataset and those recovered previously, indicating that the ndna data sufficiently overpowered the conflicting cyt b signal. the two sub-clades recovered in our analysis are congruent with the patterns detected by jenkin in mandibular morphology. the shallow-keeled clade comprises those species with bills suited to capture larger prey items and includes the greater, american and chilean flamingos. within this clade, the american and greater flamingos are often considered sub-species. the 1% divergence inferred from the dna barcode locus falls in the range of expected nearest neighbor distances among sister species, such as seen commonly in waterfowl. although more sampling is needed to confirm species status, our finding of fixed mutational differences in several nuclear loci and a 1% divergence in mtdna is consistent with treating the american and greater flamingos as separate species. the deep-keeled clade comprises those species with bills suited to capturing smaller prey items and includes the lesser, andean and james’s flamingos. these species have often been set apart from the classical genus phoenicopterus with the lesser placed within phoeniconaias based on unique characteristics of the bill and the andean and james’s placed together in phoenicoparrus based on the shared absence of the hallux. however, the morphological definition of phoeniconaias does not seem to exclude andinus and jamesi, and the hallux of minor is reduced, approximating the morphological definition of phoenicoparrus. thus, we suggest that the division of the deep-keeled clade into two genera is arbitrary and uninformative. instead, we suggest the placement of minor within phoenicoparrus along with andinus and jamesi based on shared mandibular morphology, ecology and phylogenetic relatedness. we suggest the designation of phoeniconaias as a junior synonym of phoenicoparrus based on priority and the redefinition of phoenicoparrus as all species more closely related to phoenicoparrus andinus than to phoenicopterus roseus. divergence time estimation even though the various subsets of data recovered incongruent divergence times, all agreed on a recent divergence for crown phoenicopteridae, indicating that living flamingos are among the youngest lineages within neoaves. this young age is counter to the general perception of the phoenicopteridae as among the oldest lineages of living birds based on fossil data but is congruent with the dna-dna hybridization data of sibley and ahlquist, which identified the close genetic relatedness of the flamingos and predicted an age of 5– ma for the basal divergence among the extant species. the ndna dataset recovered ages significantly younger than those of the mtdna dataset. this pattern is likely due in part to artificial signal caused by saturation of the more quickly mutating mtdna genes. though saturation is unlikely to be an issue with divergences as shallow as those of crown phoenicopteridae, the use of temporally more distant calibrations external to the clade of interest may have introduced saturation artifact. it is likely that homoplasy within the relatively longer and more quickly-evolving mtdna sequences, with more mutational sites, conflicted with the signal of the shorter and more slowly evolving ndna sequences, resulting in artificially old age estimates. thus, the true ages of the component flamingo divergences are probably somewhere between those predicted by the total dataset and by the ndna dataset alone. biogeography the two old world species are each recovered in clades with otherwise exclusively new world distributions. the recovery of the greater flamingo as derived within phoenicopterus suggests a new world origin for the shallow-keeled flamingos. the origin of phoenicoparrus is not so straightforward, as the basal divergence within the clade spans the atlantic. plio-pleistocene fossil flamingos are known only from the new world and australia and the phylogenetic affinities of these species are uncertain. no flamingos occupy australia today and it is unknown if the plio-pleistocene fossils represent an extension of the current distribution or if they are descendants of the species which occupied australia in the oligo-miocene. the new world fossil record, however, is congruent with the modern flamingo distribution and likely suggests an origin for phoenicoparrus, and crown phoenicopteridae, in the western hemisphere, followed by trans-atlantic dispersal. better phylogenetic resolution for the australian forms will further inform phoenicopterid origins. implications for the fossil record the phoenicopteriform fossil record is rich but poorly understood and revision will greatly advise divergence time estimation and paleo-biogeographical reconstruction for crown flamingos. the stem phoenicopteriform record is characterized by fragmentary material with only tentative assignment to the flamingos. while the earliest palaelodids are the earliest to be reliably placed in phoenicopteriformes, juncitarsus, elornis and agnopterus from the middle to late eocene of north america and europe may be earlier members and perhaps more appropriate fossil calibrations for mirandornithes. better understanding of the plio-pleistocene record may also provide internal calibrations for crown flamingos. finally, while we present the first information on the timing of diversification and specialization within crown flamingos, better understanding of interspecific variation is necessary to interpret the ecological significance of these results. likewise, revision of the stem phoenicopteriform fossil record is necessary to place these ages in a more meaningful context. most notably, the primitive flamingo harrisonavis croizeti from the oligo-miocene of europe displayed a mandibular morphology approximating that of modern flamingos. reassessment of this material, and new descriptions of contemporaneous taxa, will inform the state of flamingo specialization at this point and better constrain the rate of specialization within the phoenicopteridae.  chilean flamingos show no evidence of population structure across their entire range, a condition similar to other widely-distributed flamingo species. the six living species of flamingos fall into two clades reflecting differences in bill morphology and feeding ecology. we propose slight reorganization of flamingo taxonomy to indicate these differences: the greater, american and chilean flamingos comprise genus phoenicopterus while the lesser is moved to genus phoenicoparrus along with the james’s and andean flamingos. the living species diverged from each other ~ million years ago, followed by phoenicoparrus and phoenicopterus. methods and materials population genetics methods field sampling of all flamingo populations remains impractical for logistic reasons, including geographic inaccessibility and conservation protection. sampling from captive birds alone is not guaranteed to capture species-wide variation due to a persistent lack of information regarding the geographic origins of the ancestors of captive populations. to assess phylogeographic structure, we thus set out to sample museum specimens of chilean flamingos with known locality. available dna barcodes from chilean flamingos suggest the presence of at least two haplotypes in argentina but no data is available from other portions of their range. thus, we examined portions of three mitochondrial genes through amplification from historic toe pads from chilean flamingos with known localities. the samples used spanned nearly the entire range known of chilean flamingos and were obtained through several museum loans. the three loci were chosen based on known intraspecific variation in flamingos or known high rates of mutation. primers were manually generated using sequences available in genbank. in order to prevent contamination, ancient dna procedures for extraction and amplification were followed. taxonomic sampling tissue samples were obtained for each flamingo species. podilymbus podiceps and aechmophorus clarkii were included to represent the basal split within grebes. the sister taxon to remains controversial but the tropicbirds are consistently found as close relatives and were used as an outgroup. penguins and tubenoses were included as representatives of the waterbird clade. sample collection was reviewed and approved by the centre de recherche sur la biologie des populations d’oiseaux under the personal permit of alan johnson and arnaud béchet and by usda permit of mvt. genomic sampling genomic dna was extracted from tissue samples using the dneasy tissue kit and was amplified for nuclear and two mitochondrial loci from one individual per species, which comprised the primary dataset. two additional nuclear loci were also amplified to test intraspecific variation among flamingos only. see additional file for primer sequences. the nuclear loci consisted primarily of intronic regions with primers designed using an exon-primed intron spanning approach and sequences from the chicken and zebra finch genome. polymerase chain-reaction amplifications were carried out in µl reactions comprising: µl dh2o, µl 10× detergent- and mg+-free reaction buffer, µl mgcl, µl dntp, µl each primer, µl hot multitaq polymerase and µl dna template. thermocycling comprised an initial denaturation at 95°c for min followed by 35– cycles of annealing at 95°c for sec, extension at 50°c for sec and denaturation at 95°c for sec, followed by a final annealing/extension at 72°c for min. excess dntps and primers were removed from pcr product by addition of µl exonuclease, µl fastap shrimp alkaline phosphatase and µl dh and held at 37°c for min and then 80°c for min. product was sequenced offsite at macrogen usa. only the two mtdna loci were treated as coding; the mtdna were further partitioned by codon in the divergence time estimation analyses. substitution models for each locus were estimated using mega and the best available model was used in mega and mrbayes. for the phylogenetic analyses and in beast. for the divergence time estimation analyses. adamts and hmgb were only used in phylogenetic analysis testing flamingo intraspecific variation. abbreviations: rheb rheb isoform intron 3; timm17a, mitochondrial import inner membrane translocase subunit tim17a, intron 3; tcf transcription factor intron 12; rps ribosomal protein s intron 5; slc29a solute carrier family member intron 8; nfkbiz, nuclear factor of light polypeptide gene enhancer in b-cells inhibitor, zeta, intron 6; g3pdh, glyceraldehyde-3-phosphate dehydrogenase gene, intron 11; myoglobin gene, exon and intron 2; zenk, zinc finger protein, exon 2; zenk 3’utr, zinc finger protein, 3’ untranslated region; coi, cytochrome oxidase subunit i; cyt b, cytochrome b; adamts adam metallopeptidase with thrombospondin, type motif, intron 5; hmgb high mobility protein group, box intron drawings by m. mccracken. sequence data for one to two additional flamingo individuals per species for five loci were collected in the dijon, france lab. for these samples, genomic dna was extracted from blood samples for the greater, american and chilean flamingos and feather samples for the remaining three species using a standard phenol-chloroform method. pcr amplifications and product clean-up followed the same protocols as the wilmington, us lab. product was sequenced offsite at macrogen korea. sequence treatment and phylogenetic analysis sequences were aligned within sequencher and alignments were confirmed within mega using clustalw with default settings. pairwise genetic distances were calculated in mega using the default maximum composite likelihood substitution model except for coi, which was calculated using kimura 2-parameter. maximum likelihood phylogenetic analyses were carried out within mega and the loci were concatenated. unless the dataset comprised only coding genes, each dataset was treated as noncoding; optimal substitution models were found within the hierarchical bi and ml framework in mega and the analysis used a nearest-neighbor-interchange heuristic search. bootstrap support was estimated from replicates using the same settings. primary phylogenetic analyses using bayesian inference were carried out in mrbayes.. datasets were fully partitioned based on the optimal substitution models estimated within mega. each search comprised two concurrent runs of four chains each for, generations sampled every generations with the first, generations discarded as burn-in. the results of each analysis were tested for convergence of phylogenetic signal in tracer, where estimated sample size values > were treated as reliable signal. the primary analysis employed all ndna and two mtdna loci for one individual per species. all secondary analyses employed the same parameters as the primary analysis except where noted. to test for the presence of intraspecific variation, a series of single-locus analyses were carried out using two to three sequences per flamingo species for five loci. to recover individual gene trees, each locus was analyzed independently. to test the relative weight of each locus to the combined dataset, a series of analyses were run excluding a single gene at a time. to test for conflicting signal between the ndna and mtdna datasets, each was analyzed independently. to test for the effects of over-partitioning of the bayesian inference analysis, a bayesian analysis was carried out treating the total dataset as non-partitioned. a bayes factor analysis was carried out in tracer to compare the strength in phylogenetic signal of the fully partitioned and non-partitioned bi treatments. *beast. was used to test for discordance between individual gene trees and the species tree. the dataset was fully phased by individual and gene. loci were allowed unlinked substitution models, clock models and partition trees. the optimal available substitution model was used and the coding loci were partitioned by codon position. each locus was allowed a lognormal relaxed clock with uncorrelated rates. mean clock rate priors were allowed uninformative uniform distributions. a yule speciation process was assumed. no fossil node calibrations were used. the results were tested for convergence of phylogenetic signal in tracer. divergence timing analysis and fossil calibrations divergence time estimation analyses were carried out using beast.. loci were allowed unlinked substitution models and the coding loci were partitioned by codon position. a lognormal relaxed clock with uncorrelated rates was assumed and a yule speciation process was employed. the clock mean was allowed an uninformative uniform prior distribution. taxa were divided into sets representing grebes and mirandornithes; monophyly of these sets was not enforced. two nodes were fossil calibrated following the guidelines suggested by parham et al.. the divergence of extant grebes was calibrated at mya based on thiornis sociata, a grebe from the late miocene of spain. the divergence of flamingos and grebes was calibrated at mya based on adelalopus hoogbutseliensis, a stem phoenicopteriform from the earliest oligocene of belgium. calibrations were allowed lognormal distributions with offset = calibration age and standard deviation = . see table for the data included in each divergence analysis. all analyses were run twice for million generations sampled every with burn-in of million. convergence of phylogenetic signal was tested for in tracer. availability of supporting data all supporting data are included as additional files. newly obtained dna sequence reads have been uploaded to genbank. genbank accession numbers are provided in additional files and competing interests the authors declare that they have no competing interests. authors’ contributions crt carried out dna extraction, amplification, alignment, performed all phylogenetic and divergence time analysis, and drafted the manuscript. lmo and mafg carried out dna extraction, amplification, and alignment, and helped draft the manuscript. bf carried out ancient dna extraction, amplification and alignment. mvt conceived of the study, participated in its design and coordination, assisted with alignment and analysis, and helped to draft the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 specimen and locality information for the chilean flamingo individuals included in the population genetics analysis. click here for file additional file 2 pairwise genetic distances between multiple individuals of each flamingo species for five ndna loci. click here for file additional file 3 pairwise genetic distances between phoenicopterus ruber and p. roseus for ndna and two mtdna loci. click here for file additional file 4 primer sequences for each locus used in this study. click here for file additional file 5 ancient dna extraction and amplification methods. click here for file additional file 6 sample information for individuals used in the primary phylogenetic analyses. click here for file additional file 7 genbank accession numbers for sequence data used in the primary phylogenetic analyses. click here for file additional file 8 sample information for the additional flamingo individuals used in the pairwise genetic distance analyses and genbank accession numbers for associated sequence data. click here for file acknowledgements we would like to thank the many people and institutes for sample collection and/or approving sample loans: sharon birks, brooks childress and rebecca lee, james dean, jack demboski, becky desjardins and john gerwin, steve emslie, maureen flannery and jack dumbacher, blair hedges, janet hinshaw, peter houde, diana humple, irby lovette, steve cardiff, donna dittman and van remsen, mark robbins, paul sweet and joel cracraft, krystof ziskowski, sylvan heights waterfowl centre, parc des oiseaux, parc zoologique d’amnéville, salins group and tour du valat. we further like to thank frank cézilly, arnaud béchet and rémi wattier for access to samples and facilities at the université de bourgogne, dijon, france, and dan ksepka for discussion of the flamingo and grebe fossil record. for financial support, we thank the college of arts and sciences and department of biology and marine biology at uncw, the conseil régional de bourgogne and the total foundation. we thank m. mccracken for providing flamingo drawings.
modern flamingos occupy a highly specialized ecology unique among birds and represent a potentially powerful model system for informing the mechanisms by which a lineage of birds adapts and radiates. however, despite a rich fossil record and well-studied feeding morphology, molecular investigations of the evolutionary progression among modern flamingos have been limited. our divergence time estimates place flamingos among the youngest families of birds, counter to the classical notion of flamingos as among the oldest based on biogeography and the fossil record. phylogenetic analyses include likelihood and bayesian frameworks and account for potential gene tree discordance. here, using three mitochondrial markers, we present the first dna sequence-based study of population genetic variation in the widely distributed chilean flamingo and, using two mtdna and nuclear markers, recover the species tree and divergence time estimates for the six extant species of flamingos. mtdna sequences confirmed the presence of a single metapopulation represented by two minimally varying mtdna barcodes in chilean flamingos. the finding of a single metapopulation in the widespread chilean flamingo confirms similar findings in other widespread flamingo species. finally, we designate ‘phoeniconaias’ as a junior synonym of ‘phoenicoparrus’ and redefine the latter genus as containing all flamingos more closely related to phoenicoparrus andinus than phoenicopterus roseus. analyses of divergence times are fossil calibrated at the divergence of mirandornithes and the divergence of crown grebes. there is significant incongruence between the ages recovered by the mtdna and ndna datasets, likely due to mutational saturation occurring in the mtdna loci. the initial divergence among flamingos occurred at or shortly after the mio-pliocene boundary followed by quick consecutive divergences throughout the plio-pleistocene.
 modern flamingos occupy a highly specialized ecology unique among birds and represent a potentially powerful model system for informing the mechanisms by which a lineage of birds adapts and radiates. however, despite a rich fossil record and well-studied feeding morphology, molecular investigations of the evolutionary progression among modern flamingos have been limited. here, using three mitochondrial markers, we present the first dna sequence-based study of population genetic variation in the widely distributed chilean flamingo and, using two mtdna and nuclear markers, recover the species tree and divergence time estimates for the six extant species of flamingos. phylogenetic analyses include likelihood and bayesian frameworks and account for potential gene tree discordance. analyses of divergence times are fossil calibrated at the divergence of mirandornithes and the divergence of crown grebes.  mtdna sequences confirmed the presence of a single metapopulation represented by two minimally varying mtdna barcodes in chilean flamingos. likelihood and bayesian methods recovered identical phylogenies with flamingos falling into shallow-keeled and deep-keeled sub-clades. the initial divergence among flamingos occurred at or shortly after the mio-pliocene boundary followed by quick consecutive divergences throughout the plio-pleistocene. there is significant incongruence between the ages recovered by the mtdna and ndna datasets, likely due to mutational saturation occurring in the mtdna loci.  the finding of a single metapopulation in the widespread chilean flamingo confirms similar findings in other widespread flamingo species. the robust species phylogeny is congruent with previous classifications of flamingos based on feeding morphology. modern phoenicopterids likely originated in the new world with each sub-clade dispersing across the atlantic at least once. our divergence time estimates place flamingos among the youngest families of birds, counter to the classical notion of flamingos as among the oldest based on biogeography and the fossil record. finally, we designate ‘phoeniconaias’ as a junior synonym of ‘phoenicoparrus’ and redefine the latter genus as containing all flamingos more closely related to phoenicoparrus andinus than phoenicopterus roseus. flamingophylogenydivergence timebiogeographyfossilbillfilter feeding 