Dna barcoding, population genetics, and phylogenetics of the illegally hunted Philippine Duck Anas luzonica (Aves: Anseriformes: Anatidae)

 

 

 

 

Article

INTRODUCTION

Bird taxonomy during the greater part of the 20th century was the subject of much debate. Although several morphology-based phylogenetic studies investigated the evolutionary relationships of taxa within class Aves, findings were constantly challenged (Livezey & Zusi 2007). The application of molecular biology techniques to avian systematics finally gave way to the wide acceptance of hypotheses regarding divergences at the base of the avian tree (Meyer & Zardoya 2003; Hackett et al. 2008). Higher-order relationships remain largely unresolved, however (Livezey & Zusi 2007).

It was only recently that a robust genus-level phylogeny for order Anseriformes (screamers and waterfowl) was generated. The robust tree was constructed based on the analysis of sequences of the complete mitochondrial DNA control region (Donne-Gousse et al. 2002). Resolving relationships deeper within this order, especially those among the dabbling ducks (Anatidae: Anas) requires more complicated methods. This is because allele sharing likely occurs among the dabbling ducks due to the relatively recent emergence of this lineage and frequent hybridization among its members (Johnson & Sorenson 1999; Lavretsky et al. 2014). To address these complications, Lavretsky et al. (2014) generated species trees based on 20 nuclear loci using multilocus coalescent methods. The nuclear topologies were compared to mitochondrial and morphometric topologies. Generally, the major clades were composed of the same taxa across all topologies, resolving the relationships among the said taxa. The only taxa whose placement within the different trees varied were the Hawaiian Duck, the Indian and Chinese Spot-billed Ducks, and the Philippine Duck (Lavretsky et al. 2014).

The Philippine Duck Anas luzonica is the only endemic duck found in the Philippines (Kennedy et al. 2000). It is listed as ‘Vulnerable’ by the International Union for Conservation of Nature (IUCN) mainly due to rapid population declines caused by continuous hunting pressures. In 2002, the total population of Anas luzonica was estimated to be fewer than 10,000 individuals (BirdLife International 2012). Despite being protected by law, Philippine Ducks continue to be hunted for sport and for their meat (Galang 2014). Since duck meat can easily be missed by authorities when labeled by illegal hunters or smugglers as domestic meat, there is a need for a method of identification which does not rely on the morphological characters of the Philippine Duck. DNA barcoding may fill this need. Determining the DNA sequence from a standard part of the genome of the unknown specimen and comparing it to other existing ‘barcodes’ provides a means of species identification (Hajibabaei et al. 2007). The mitochondrial cytochrome c oxidase subunit I (COI) gene is the most commonly used marker for barcoding animals, since it possesses a relatively strong phylogenetic signal and its primers are effective in recovering the gene in almost all animal phyla (Hebert et al. 2003).

Using muscle tissue samples obtained from illegally hunted and confiscated Philippine Ducks from Pantabangan, Nueva Ecija, the present study aimed to generate a COI barcode profile for Anas luzonica. The study also aimed to examine the population genetics of Philippine Ducks using mitochondrial DNA control region and COI sequences. Finally, given the lack of certainty regarding the relationship of Anas luzonica with other dabbling ducks, the study aimed to generate a phylogenetic tree for Anas spp. based on the COI gene.

 

 

MATERIALS AND METHODS

 

 

Sampling

Muscle and liver tissue samples were obtained from 46 individuals of Philippine Duck from the town of Pantabangan, Nueva Ecija (15.81670N & 121.15000E). The ducks were confiscated by staff of the Department of Environment and Natural Resources Region III Office on 12 February 2014 from illegal hunters who targeted the IUCN Red Listed Vulnerable species.

 

DNA extraction

The modified NaOH-lysis method (Fontanilla 2010) was employed to extract DNA from muscle tissue samples of 46 individuals of Anas luzonica. The tissues were sliced using sterile glass slides, then transferred to 1.5-ml microcentrifuge tubes. The slices were then ground in 400μl of 0.1 N NaOH using glass beads and a plastic pestle. When pieces of tissue were no longer visible, the samples were heated at 95­­–100°C for 20 min in a dry bath. Afterwards, 200μl of sterile distilled water and 600μl of chloroform-isoamyl alchohol (24:1) were added. The tubes were inverted 100 times, then centrifuged at 13,000 rpm for 10 min. The upper phase was transferred to new tubes, and an equal volume of isopropanol was added. The tubes were inverted several times before being stored at -80°C for one hour. The samples were spun at 13,000 rpm for 15 min, after which the isopropanol was discarded, leaving the pellets in the tubes. These were washed with 500μl of cold 70% ethanol, then spun at 13,000 rpm for five minutes. The ethanol was decanted out. To ensure that the pellets were dry, the tubes were heated at 450C for 15 min at most. The pellets were then suspended in 150μl of ultrapure water. The resulting DNA extracts were stored at -200C until use.

 

PCR amplification

The COI gene was amplified via the polymerase chain reaction (PCR) using primers BirdF1 (5’-TTCTCCAACCACAAAGACATTGGCAC-3’) and BirdR1 (5’-ACGTGGGAGATAATTCCAAATCCTG-3’) (Hebert et al. 2004). Fifty μl reactions were prepared in 0.2-ml PCR tubes. The components of each reaction were as follows: five μl of 10X PCR buffer, one μl of 1.25mM dNTPs, 10 μl of QIAGENTM Q solution, two μl of 25mM MgCl2, 2.5μl of each primer (10μM), 22.75 μl of ultrapure water, 0.25μl of RocheTM 5 units/μl Taq polymerase, and four μl of the DNA template (0.5-142.9 ng/µl). PCR reactions were carried out in a Labnet MultiGeneTM thermocycler under the following conditions: 60 sec of initial denaturation at 940C, followed by six cycles of 60 sec at 940C, 90 sec at 450C, and 90 sec at 72°C, then 35 cycles of 60 sec at 940C, 90 sec at 550C and 90 sec at 720C, and a final extension of five minutes at 720C (Kerr et al. 2007).

To study the population genetic structure of the Anas luzonica population, the 5’ half of the mitochondrial DNA control region (bp 78–774) was also amplified (Desjardins & Morais 1990). The primers used to amplify this region are L78 (5’- GTTATTTGGTTATGCATATCGTG-3’) and H774 (5’-CCATATACGCCAACCGTCTC-3’) (Sorenson & Fleischer 1996; Sorenson et al. 1999). Fifty μl reactions, composed of five μl of 10X PCR buffer, one μl of 1.25mM dNTPs, 2.5μl of each primer (10 μM), 34.75μl of ultrapure water, 0.25μl of RocheTM5 units/μl Taq polymerase, and four μl of the DNA template (0.5–142.9 ng/µl), were prepared. Thermocycling was conducted under the following conditions: seven minutes of initial denaturation at 940C, 38 cycles of 20 sec at 940C, 20 sec at 520C, 60 sec at 720C, and a final extension of seven minutes at 720C (McCracken et al. 2001).

 

PCR product purification and sequencing

To ensure the absence of contamination, the PCR products were electrophoresed in 1% agarose gel at 100 volts for 30 min, with 0.5X TBE (Tris base, boric acid, and ethylene diamine tetra-acetic acid) as the running buffer. The in-gel staining method, using 1% ethidium bromide, was used for visualization under ultraviolet illumination. The resulting bands were excised and placed in 2-ml microcentrifuge tubes. These were stored at -200C until further processing.

The QIAGENTM QIAquick Gel Extraction Kit (Germany) was used following the manufacturer’s protocol to purify the PCR products from the excised gels. The products were then sent to 1st Base in Singapore for Sanger dideoxy sequencing using the same forward and reverse primers used in PCR amplification to sequence both strands of each region.

Sequence assembly and alignment

The Anas luzonica DNA sequences were assembled using the STADEN package (Staden et al. 2000). The COI gene sequences of other Anas species along with those of outgroup taxa in the subfamily Anatidae (Aix galericulata, Aix sponsa, Aythya affinis, Aythya americana, Aythya collaris, Aythya marila, Aythya valisineria, Callonetta leucophrys, Chenonetta jubata, Sarkidiornis melanotos, Tadorna ferruginea, and Tadorna tadorna) were obtained from GenBank (www.ncbi.nlm.nih.gov/genbank/), available as of 10 June 2014 (see Table 1). A total of 264 COI sequences, including the 46 Anas luzonica COI sequences from this study, were aligned using Clustal W in the BioEdit Sequence Alignment Editor (Hall 1999) and adjusted manually by eye. The same software was used to align 45 Anas luzonica control region sequences. Sequences were trimmed to limit nucleotide positions to only those common to all sequences. The COI sequences were trimmed to a length of 610 bp, while the control region sequences were trimmed to 622 bp.

All A. luzonica sequences were submitted to GenBank with the following accession numbers: KT151676 - KT151721 for COI and KT151722 - KT151766 for the mitochondrial control region.

Sequence analysis

The optimal model of DNA substitution was determined using the Bayesian Inference Criterion (BIC) method under the jModelTest ver. 2.1.7 program (Posada & Crandall 2001). Using PAUP* version 4.0b10 (Swofford 2002), pairwise genetic distances based on the optimal model were estimated. The maximum likelihood method (Cavalli-Sforza & Edwards 1967) was applied to construct a COI gene tree using MEGA version 7 (Kumar et al. 2016). Bootstrap resampling (Felsenstein 1985) was performed to quantify the support for the branches of the generated tree.

Analysis of population genetic structure

Both COI and control region sequences were used for population genetic analysis. DnaSP version 5.10 (Librado & Rozas 2009) was used to estimate the haplotype diversity (h) and nucleotide diversity (p) of the Anas luzonica population based on Nei & Li’s (1979) formulas.

Haplotype diversity (h) is commonly used as a measure of population variation. It is the probability of two haplotypes, randomly drawn from the population, being different (Nei 1987). Nucleotide diversity (p) is the average number of nucleotide differences per site between two DNA sequences selected at random.

RESULTS AND DISCUSSION

 

 

This study successfully generated 46 COI sequences from Anas luzonica, with three distinct haplotypes found. The division of the 46 COI sequences into the three haplotypes are as follows: 27 for haplotype 1, 18 for haplotype 2, and one for haplotype 3. Haplotypes 2 and 3 differ from haplotype 1 in only one nucleotide position. Haplotype 2 differs from haplotype 3 in two positions. The nucleotide positions where the haplotypes vary are shown in Table 2.

For the genetic distance and neighbor joining analyses, a total of 26 Anas species (including A. luzonica) were evaluated based on the availability of their COI sequences in Genbank. In the species considered, the number of sequences analyzed and the number of unique sequences varied (Table 3). Overall, 80 unique sequences were identified from the 264 sequences used.

Genetic divergences were determined using the TrN+Γ+I model of DNA evolution, which was determined to be the optimal model based on the BIC criterion in jModelTest. Genetic divergences, measured as the average TrN+Γ+I distance, increased as the taxonomic level of the comparisons became more inclusive (Table 4). The mean TrN+Γ+I distance between genera (17.33%) was more than double that between species (8.24%), while the latter was nearly 55 times higher than the mean conspecific TrN+Γ+I distance value (0.15%).

The effectivity of barcoding in species delineation is dependent upon the presence of a barcode gap, i.e. the separation between the range of intra- and interspecific variations (Meyer & Paulay 2005). In the COI sequences analyzed, the highest intraspecific TrN+Γ+I distance (7.24%) exceeded the smallest interspecific TrN+Γ+I distance (0.16%), indicating the absence of the gap. Figure 1 illustrates the degree of overlap between intra- and interspecies variation. Majority (99.73%) of variation within species was found to occur within the 1.3% distance based on the TrN+Γ+I model, which is less than the value derived by Hebert et al. (2004) for birds (2.7%) based on the less realistic K2P model. Some (10.44%) of the comparisons between species however also fell within this distance, with 2.33% of interspecies comparisons showing very little genetic divergence.

The utility of DNA barcoding for distinguishing Anas luzonica from other Anas species was supported by the maximum likelihood tree constructed based on COI gene sequences (Fig. 2). All three Anas luzonica haplotypes clustered into a monophyletic group supported by a bootstrap value of 96%. This monophyly may be considered significant since Anas luzonica is known to hybridize with Anas platyrhynchos and Anas superciliosa (Pacific Black Duck), both of which also belong to the mallard clade (McCarthy 2006). Sampling of Anas luzonica from geographically widespread locations, especially areas of sympatry with the Mallard and Pacific Black Duck, may provide further insight into the utility of DNA barcoding for Anas luzonica species delineation.

Based on the COI sequences and the NJ tree, this study noted Anas species with identical haplotypes such as A. cyanoptera and A. discors; A. americana, A. flavirostris and A. sibilatrix; A. platyrhynchos and A. poecilorhyncha; and A. platyrhynchos, A. poecilorhyncha and A. rubripes. Sharing of haplotypes may explain why some species lump together in the COI maximum likelihood tree (Fig. 2). Haplotype sharing between Anas species is not an uncommon occurrence. Such has been observed between the mitochondrial control region sequences of A. penelope and A. americana, and those of A. penelope and A. strepera. Each of these pairings is known to hybridize and produce viable offspring. Haplotype sharing between these species has been attributed to the introgression of mitochondrial DNA into previously genomically distinct species (Peters et al. 2005). Interspecies hybridization may also be the cause of the lack of genetic divergence between some species in the present study. Cases of hybridization in the wild have been reported for Anas cyanoptera (Cinnamon Teal) and Anas discors (Blue-winged Teal), which exhibit paraphyly with respect to each other within the clade of blue-winged ducks (Cockrum 1952; Kessler & Avise 1984; Johnson & Sorenson 1999); however, the paraphyly of these two species may have also been a result of a relatively recent divergence, meaning A. cyanoptera and A. discors are young lineages and ancestral haplotypes are still retained. The cause of low divergence and paraphyly in A. cyanoptera and A. discors warrants further study (Johnson & Sorenson 1999).

Following its widespread geographic distribution, the Mallard (A. platyrhynchos) has been well known to hybridize and produce 100% viable offspring with other species, including the American Black Duck (A. rubripes) (Avise et al. 1990). Before the European settlement in North America, A. platyrhynchos and A. rubripes populations were allopatrically separated. The latter was considered an isolated offshoot of the former. Habitat alteration and game-farm mallard releases during the 20th century resulted in sympatry between the species. Increased secondary contact has led to extensive hybridization between A. platyrhynchos and A. rubripes (Mank et al. 2004). Mank et al. (2004) speculated that hybridization brought about the significant reduction in genetic differentiation for Mallards and Black Ducks. A later study which examined the effect of the exclusion/inclusion of introgressing taxa on relationships among the Mallards and allies revealed that at least for New World taxa (which includes A. rubripes), hybridization plays a significant role in confounding phylogenetic relationships (Lavretsky et al. 2014).

In the present study, the other species which shared COI haplotypes with the Mallard was the Indian Spot-billed Duck A. poecilorhyncha. Although there is a lack of documentation for fertile hybrids formed between A. platyrhynchos and A. poecilorhyncha individuals, A. platyrhynchos is known to hybridize with the Chinese Spot-billed Duck A. zonorhyncha, the Indian Spot-billed Duck’s eastern counterpart (Kulikova et al. 2004). Upon Random Amplified Polymorphic DNA analysis, a study by Kulikova et al. (2003) reported low genetic differentiation between A. platyrhynchos and A. poecilorhyncha. The geographic range of Spot-billed Duck described in the study, which covered northern and western China, as well as Korea and Japan, indicates that the species is likely to be A. zonorhyncha rather than A. poecilorhyncha (Kulikova et al. 2003).

Cross accounts of A. americana American Wigeon with A. sibilatrix Chiloé Wigeon have been documented; however, these hybridization events have occurred in captivity as A. americana is primarily restricted to North and Central America, while A. sibilatrix mainly inhabits the southern cone of South America (McCarthy 2006; Kerr et al. 2009). Despite differences in plumage characteristics and behavioral displays, COI sequence overlaps between the two species have been recorded (Kerr et al. 2009). Peters et al. (2005) explained that the rapid accumulation of phenotypic differences is a probable result of differing selective pressures between the northern and southern hemispheres (Kerr et al. 2009).

Haplotype sharing of A. flavirostris Speckled Teal with both A. americana and A. sibilatrix is difficult to explain because the taxa belong to different clades (Fig. 2). Moreover, there are no documented hybridization events between the species.

With the exception of one haplotype of A. flavirostris, all species analyzed grouped in accordance to their respective clades (Johnson & Sorenson 1999). Monophyly was exhibited by mallards, pintails, Baikal teals, and blue-winged ducks. Some of the phylogenetic relationships between clades were similar to those elucidated by Johnson & Sorenson (1999) using cytochrome-b and NADH dehydrogenase subunit 2 (ND2) mitochondrial DNA sequences, with some notable exceptions. In the current study, for example, the silver teals were not monophyletic, with A. quercedula falling at the base of the blue-winged ducks with 82% bootstrap support.

Results obtained by Johnson & Sorenson (1999) and the present study agree that for the wigeon/pintail/teal/mallard clade, the lineage of wigeons branched out early (100% bootstrap support in the cyt b and ND2 consensus tree and 52% in the current COI maximum likelihood tree). The sister clade of the mallards differed however; for the cyt b and ND2 consensus tree, this was the green-winged teals, while in the current COI maximum likelihood tree, it was the grey teal (90% bootstrap support).

The discrepancies between the COI NJ tree generated in the study and trees constructed based on other genes support the argument made by Lavretsky et al. (2014) that multiple genes need to be considered in resolving relationships among dabbling ducks, since extensive allele sharing occurs within this taxon.

Haplotype and nucleotide diversities for COI and control region sequences of the A. luzonica population are listed in Table 5. Diversities for the COI gene are based on 46 sequences while those for the control region are based on 45 sequences. One control region sequence was poor in quality and was deemed unusable for analysis.

Haplotype diversities of both gene sequences are high (h>0.5), while the nucleotide diversities are low (p<0.005). According to Grant & Bowen (1998), stable populations with large population sizes are characterized by high haplotype and nucleotide diversities (h>0.5, p>0.005). High h and low p may be a consequence of the rapid expansion of a population after the occurrence of a bottleneck event in which the effective population size is significantly reduced. The high haplotype diversity is explained by the enhancement in retention and accumulation of new mutations when population growth following the bottleneck event is rapid (Grant& Bowen 1998).

Following the interpretation proposed by Grant & Bowen (1998) on high haplotype diversity and low nucleotide diversity, the A. luzonica population investigated appears to be undergoing a population bottleneck. The bottleneck effect, the drastic reduction in size of a population, may occur when a considerable number of individuals are killed by a natural calamity (e.g., fire, earthquake, flood, etc.) or by human activities (Campbell et al. 2008). Human activities are the likely causes of the bottleneck effects experienced by the population under study. Utilizing other markers, particularly those from the nuclear genome, may confirm the bottleneck phenomenon in the Philippine Duck.

Major threats to populations of the Philippine Duck include hunting and trapping. The excessive practice of these activities, which first became evident in the 1960s, persists to this day, as evidenced by the fact that the tissue samples used in the present study were confiscated from illegal hunters (BirdLife International 2012). Populations are also impacted indirectly when their wetland habitats are disturbed. For example, draining wetland ecosystems for irrigation and agricultural reclamation leaves little habitat for waterfowl (Department of Environment and Natural Resources & United Nations Environment Programme 1997; Crosby 2003). The use of wetlands for aquaculture, i.e., the conversion of these ecosystems to shrimp and fish ponds, has caused near-extirpations of local populations of the Philippine Duck. The destruction of mangrove areas, which Philippine Ducks also inhabit, is associated with reduced population sizes as well. Lastly, runoff rich in pesticides from rice fields is also considered a major threat to populations of the Philippine Duck (BirdLife International 2012).

Although its range covers at least 30 localities in the Philippines extending from Luzon to Mindanao, the Philippine Duck is legally protected only in five localities, namely Mounts Iglit-Baco and Naujan Lake National Parks in Mindoro, Aurora Memorial National Park, Bataan Natural Park and Subic Bay Forest Reserve, and Olango Island (BirdLife International 2012). Technically, the Philippine Duck is protected by the Wildlife Resources Conservation and Protection Act or Republic Act No. 9147 (Department of Environment and Natural Resources 2015). However, since the passage of the Act on July 2001, the population trend of the Philippine Duck has continued to decrease. In 2004, the Asian Waterbird Census counted 4,632 individuals. The following year, 4,428 were counted (BirdLife International 2012).

Considering the number of threats A. luzonica faces, and the ineffective protection it has been given, it is likely that other populations of the Philippine Duck are experiencing bottleneck effects. Populations that have undergone bottleneck events may have very low levels of genetic variation for long periods of time as a consequence of reduced effective population size and inbreeding (Campbell et al. 2008). Low genetic variation does not bode well for the fitness of a population as it reduces fertility and increases embryonic mortality (Keller & Waller 2002).

The high haplotype diversities calculated in this study suggest that low genetic variation is not an issue for the A. luzonica population assessed. However, this population may not be representative of all populations of A. luzonica in the Philippines. Extensive sampling from several localities must be conducted to ascertain whether other Philippine Duck populations are at risk for inbreeding depression (Keller & Waller 2002).

 

REFERENCES

 

 

 

Avise, J.C., D.C. Ankney & W.S. Nelson (1990). Mitochondrial gene trees and the evolutionary relationship of Mallard and Black Ducks. Evolution 44: 1109–1119; http://doi.org/10.2307/2409570

BirdLife International (2012). Anas luzonica. In: IUCN 2014. 2014 IUCN Red List of Threatened Species. Downloaded on 2 July 2014.

BirdLife International (2016). Anas luzonica. The IUCN Red List of Threatened Species 2016: e.T22680214A92849560. Downloaded on 13 May 2017.

Campbell, N.A., J.B. Reece, L.A. Urry, M.L. Cain, S.A. Wasserman, P.V. Minorsky & R.B. Jackson (2008). Biology, 8th Edition. Pearson Benjamin Cummings, San Francisco, 1267pp. 

Cavalli-Sforza, L.L. & A.W.F. Edwards (1967). Phylogenetic analysis: models of estimation procedures. American Journal of Human Genetics 19: 233–257.

Cockrum, E.L. (1952). A check-list and bibliography of hybrid birds in North America, North of Mexico. The Wilson Bulletin 64: 140–159.

Crosby, M.J. (2003). Saving Asia’s Threatened Birds: A Guide for Government and Civil Society. BirdLife International, Cambridge, 256pp.

Department of Environment and Natural Resources (2015). The Philippine Duck. Retrieved from http://www.denr.gov.ph/index.php/component/content/article/621.html on 3 April 2015.

Department of Environment and Natural Resources and United Nations Environment Programme (1997). Philippine Biodiversity: An Assessment and Action Plan. Bookmark Inc., Makati City, 298pp.

Desjardins, P. & R. Morais (1990). Sequence and gene organization of the chicken mitochondrial genome: a novel gene order in higher vertebrates. Journal of Molecular Biology 212: 599–634; http://doi.org/10.1016/0022-2836(90)90225-B

Donne-Gousse, C., V. Laudet & C. Hanni (2002). A molecular phylogeny of anseriformes based on mitochondrial DNA analysis. Molecular Phylogenetics and Evolution 23: 339–356; http://doi.org/10.1016/S1055-7903(02)00019-2

Felsenstein, J. (1985). Confidence limits on phylogenies: an approach using the bootstrap. Evolution 39: 783–791; http://doi.org/10.2307/2408678

Fontanilla, I.K.C. (2010). Achatina (Lissachatatina) fulica Bowdich: its molecular phylogeny, genetic variation in global population, and its possible role in the spread of the rat lungworm Angiostrongylus cantonensis (Chen). PhD Thesis. University of Nottingham, 617pp.

Galang, A. (2014). “7 Nabbed for Hunting Wild Ducks in Nueva Ecija.” Inquirer News. Philippine Daily Inquirer, 14 Feb. 2014. Web. 03 Aug. 2014. http://newsinfo.inquirer.net/577833/7-nabbed-for-hunting-wild-ducks-in-nueva-ecija

Grant, W.S. & B.W. Bowen (1998). Shallow population histories in deep evolutionary lineages of marine fishes: insights from sardines and anchovies and lessons for conservation. Journal of Heredity 89: 415–426; http://doi.org/10.1093/jhered/89.5.415

Hackett, S.J., R.T. Kimball, S. Reddy, R.C.K. Bowie, E.L. Braun, M.J. Braun, J.L. Chojnowski, W.A. Cox, K. Han, J. Harshman, C.J. Huddleston, B.D. Marks, K.J. Miglia, W.S. Moore, F.H. Sheldon, D.W. Steadman, C.C. Witt & T. Yuri (2008). A phylogenomic study of birds reveals their evolutionary history. Science 320: 1763–1768; http://doi.org/10.1126/science.1157704

Hajibabaei, M., G.A.C. Singer, P.N.D. Hebert & D.A. Hickey (2007). DNA barcoding: how it complements taxonomy, molecular phylogenetics and population genetics. Trends in Genetics 23: 167–172; http://doi.org/10.1016/j.tig.2007.02.001

Hall, T.A. (1999). BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symposium 41: 95–98.

Hebert, P.D.N., A. Cywinska, S.L. Ball & J.R. deWaard (2003). Biological identifications through DNA barcodes. Proceedings of the Royal Society 270: 313–321; http://doi.org/10.1098/rspb.2002.2218

Hebert, P.D.N, M.Y. Stoeckle, T.S. Zemlak & C.M. Francis (2004).Identification of birds through DNA barcodes. Public Library of Science Biology 2: e312; http://doi.org/10.1371/journal.pbio.0020312

Johnson, K.P. & M.D. Sorenson (1999). Phylogeny and biogeography of Dabbling Ducks (genus: Anas): a comparison of molecular and morphological evidence. The Auk 116: 792–805.

Keller, L.F. & D.M. Waller (2002). Inbreeding effects in wild populations. Trends in Ecology & Evolution 17: 230–241; http://doi.org/10.5167/uzh-53672

Kennedy, R.S., P.C. Gonzales, E.C. Dickinson, H.C. Miranda, Jr. & T.H. Fisher (2000). A Guide to the Birds of the Philippines. Oxford University Press, New York, 369pp. 

Kerr, K.C.R., D.A. Lijtmaer, A.S. Barreira, P.D.N. Hebert & P.L. Tubaro (2009). Probing evolutionary patterns in neotropical birds through DNA barcodes. Public Library of ScienceOne 4: e4379; http://doi.org/10.1371/journal.pone.0004379

Kerr, K.C.R., M.Y. Stoeckle, C.J. Dove, L.A. Weigt, C.M. Francis & P.D.N. Hebert (2007). Comprehensive DNA barcode coverage of North American birds. Molecular Ecology Notes 7: 535–543; http://doi.org/10.1111/j.1471-8286.2007.01670.x

Kessler, L.G. & J.C. Avise (1984). Systematic relationships among waterfowl (Anatini) inferred from restriction endonuclease analysis of mitochondrial DNA. Systematic Zoology 33: 370–380.

Kulikova, I.V., G.N. Chelomina & Y.N. Zhuravlev (2003). Low genetic differentiation of and close evolutionary relationships between Anas platyrhynchos and Anas poecilorhyncha: RAPD-PCR evidence. Russian Journal of Genetics 39: 1143–1151; http://doi.org/10.1023/A:1026174910872

Kulikova, I.V., Y.N. Zhuravlev & K.G. McCracken (2004). Asymmetric hybridization and sex-biased gene flow between Eastern Spot-billed Ducks (Anas zonorhyncha) and Mallards (A. platyrhynchos) in the Russian far east. The Auk 121: 930–949; http://doi.org/10.1642/0004-8038(2004)121[0930:AHASGF]2.0.CO;2

Kumar, S., G. Stecher & K. Tamura (2016). MEGA7: Molecular Evolutionary Genetics Analysis version 7.0 for bigger datasets. Molecular Biology and Evolution 33: 1870–1874; https://doi.org/10.1093/molbev/msw054

Lavretsky, P., K.G. McCracken & J.L. Peters (2014). Phylogenetics of a recent radiation in the Mallards and allies (Aves: Anas): inferences from a genomic transect and the multispecies coalescent. Molecular Phylogenetics and Evolution 70: 402–411; http://doi.org/10.1016/j.ympev.2013.08.008

Librado, P. & J. Rozas (2009). DnaSP v5: A software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25: 1451–1452; http://doi.org/10.1093/bioinformatics/btp187

Livezey, B.C., R.L. Zusi (2007). Higher-order phylogeny of modern birds (Theropoda, Aves: Neornithes) based on comparative anatomy. II. Analysis and discussion. Zoological Journal of the Linnean Society 149: 1–95; http://doi.org/10.1111/j.1096-3642.2006.00293.x

Mank, J.E., J.E. Carlson & M.C. Brittingham (2004). A century of hybridization: decreasing genetic distance between American Black Ducks and Mallards. Conservation Genetics 5: 395–403; http://doi.org/10.1023/B:COGE.0000031139.55389.b1

McCarthy, E.M. (2006). Handbook of Avian Hybrids of the World. Oxford University Press, New York, 601pp.

McCracken, K.G., W.P. Johnson & F.H. Sheldon (2001). Molecular population genetics, phylogeography, and conservation biology of the Mottled Duck (Anas fulvigula). Conservation Genetics 2: 87–102; http://doi.org/10.1023/A:1011858312115

Meyer, A. & R. Zardoya (2003). Recent advances in the (molecular) phylogeny of vertebrates. Annual Review of Ecology, Evolution, and Systematics 34: 311–338; http://doi.org/10.1146/annurev.ecolsys.34.011802.132351

Meyer, C.P. & G. Paulay (2005). DNA barcoding: error rate based on comprehensive sampling. Public Library of Science Biology 3: e422; http://doi.org/10.1371/journal.pbio.0030422

Nei, M. (1987). Molecular Evolutionary Genetics. Columbia University Press, New York, 512pp.

Nei, M. & W. Li (1979). Mathematical model for studying genetic variation in terms of restriction nucleases. Proceedings of the National Academy of Sciences 76: 5269–5273.

Peters, J.L., K.G. McCracken, Y.N. Zhuravlev, Y. Lu, R.E. Wilson, K.P. Johnson & K.E. Omland (2005). Phylogenetics of wigeons and allies (Anatidae: Anas): the importance of sampling multiple loci and multiple individuals. Molecular Phylogenetics and Evolution 35: 209–224; http://doi.org/10.1016/j.ympev.2004.12.017

Posada, D. & K.A. Crandall (2001). Selecting the best-fit model of nucleotide substitution. Systematic Biology 50(4): 580-601; https://doi.org/10.1080/10635150118469

Sorenson, M.D., A. Cooper, E.E. Paxinos, T.W. Quinn, H.F. James, S.L. Olson & R.C. Fleischer (1999). Relationships of the extinct moanalos, flightless Hawaiian waterfowl, based on ancient DNA. Proceedings of the Royal Society of London 266: 2187–2193; http://doi.org/10.1098/rspb.1999.0907

Sorenson, M.D. & R.C. Fleischer (1996). Multiple independent transpositions of mitochondrial DNA control region sequences to the nucleus. Proceedings of the National Academy of Science USA 93: 15239–15243.

Staden, R., K. Beal & J. Bonfield (2000). The Staden package. Methods in Molecular Biology 132: 115–130.

Swofford, D.L. (2002). PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). Version 4. Sinauer Associates, Sunderland, Massachusetts, 142pp.