Incomplete lineage sorting and introgression in the diversification of Chinese spot-billed ducks and mallards.

Incomplete lineage sorting and introgression are 2 major and nonexclusive causes of species-level non-monophyly. Distinguishing between these 2 processes is notoriously difficult because they can generate similar genetic signatures. Previous studies have suggested that 2 closely related duck species, the Chinese spot-billed duck Anas zonorhyncha and the mallard A. platyrhynchos were polyphyletically intermixed. Here, we utilized a wide geographical sampling, multilocus data and a coalescent-based model to revisit this system. Our study confirms the finding that Chinese spot-billed ducks and Mallards are not monophyletic. There was no apparent interspecific differentiation across loci except those at the mitochondrial DNA (mtDNA) control region and the Z chromosome (CHD1Z). Based on an isolation-with-migration model and the geographical distribution of lineages, we suggest that both introgression and incomplete lineage sorting might contribute to the observed non-monophyly of the 2 closely related duck species. The mtDNA introgression was asymmetric, with high gene flow from Chinese spot-billed ducks to Mallards and negligible gene flow in the opposite direction. Given that the 2 duck species are phenotypically distinctive but weakly genetically differentiated, future work based on genome-scale data is necessary to uncover genomic regions that are involved in divergence, and this work may provide further insights into the evolutionary histories of the 2 species and other waterfowls.

Incomplete lineage sorting and introgression are 2 major and nonexclusive causes of species-level non-monophyly (Funk and Omland 2003; McKay and Zink 2010). Introgression can occur when hybrids are viable and fertile and backcross with at least one of the parental species (Anderson 1953). Introgression can yield species-level non-monophyly by introducing alleles across species boundaries (Funk and Omland 2003). Alternatively, incomplete lineage sorting can result in species-level non-monophyly if species divergence was too recent for ancestral polymorphisms to have sorted into reciprocally monophyletic lineages (Funk and Omland 2003; McKay and Zink 2010). Distinguishing between the 2 processes is difficult because they leave similar genetic signatures (Buckley et al. 2006; Peters et al. 2007; Wang et al. 2014).

Hybridization is common in birds, with nearly 1 in 10 species known to hybridize (Grant and Grant 1992; McCarthy 2006). The mallard complex (Anas platyrhynchos and allies), containing 11–13 recently radiated taxa (depending on taxonomic authority; Palmer 1976; Johnson and Sorenson 1999; Lavretsky et al. 2014b), is one of the best known examples with respect to the frequently documented hybridizations between Mallard and its closely related species (Johnson and Sorenson 1999; Tubaro and Lijtmaer 2002; Tracey et al. 2008; Fowler et al. 2009; Lavretsky et al. 2014b). Due to the frequent hybridization and incipient speciation, interspecific relationships within the Mallard complex are difficult to resolve (Johnson and Sorenson 1999; Gonzalez et al. 2009; Lavretsky et al. 2014b). Moreover, compared with other taxa in North America and Australasia (Avise et al. 1990; Rhymer et al. 1994; Williams et al. 2005; Guay et al. 2014; Lavretsky et al. 2014a, 2015a), the interspecific relationship between Mallard and one of its Asian relatives, the Chinese spot-billed duck (A. zonorhyncha), is poorly studied (Kulikova et al. 2003).

With a Holarctic distribution, the dichromatic Mallard’s Palearctic distribution extends from the British Isles through Europe and Russia to northern Japan and northeastern China (Carboneras and Kirwan 2018a; Figure 1). In contrast, the monochromatic Chinese spot-billed ducks have a smaller breeding range, which covers southeastern Russia, eastern Mongolia, Sakhalin, Korea, Japan, and eastern China (del Hoyo et al. 2018). Chinese spot-billed ducks are parapatrically distributed with Mallards, with distributions partially overlapping in southeastern Russia, eastern Mongolia, northeastern China, southern Sakhalin, and northern Japan (Kulikova et al. 2004).

Figure 1.

Sampling sites and distribution of mtDNA lineages (group A and group SB). Sequences of Chinese spot-billed ducks from Primorye (PR) and sequences of Old Word Mallards were downloaded from GenBank. The yellow area represents the distribution of Chinese spot-billed ducks, the purple area represents the distribution of Mallards, and the brown area represents the area of overlap between the 2 species. Circles with black borders and white borders show sampling sites of Chinese spot-billed ducks and Mallards, respectively. The circles fill colors correspond to mtDNA origin, with blue for group A haplotypes and red for group SB haplotypes. The circle size corresponds to sample size. Triangles show those sampling sites of Chinese spot-billed ducks from which samples failed to amplify mtDNA.

Previous studies suggested that Chinese spot-billed ducks and Mallards were polyphyletically intermixed (Kulikova et al. 2003, 2004; Jin et al. 2014; Lavretsky et al. 2014b) and that introgression might be the principle cause of this polyphyly (Kulikova et al. 2004). However, these results might be biased due to some constraints. First, these previous studies are limited in sample sizes and geographical coverage. Jin et al. (2014) and Lavretsky et al. (2014b) used fewer than 10 samples from each species. Kulikova et al. (2004) used a greater number of samples, but their sample locations were restricted to Primorye, Russia. Second, a limited number of nuclear loci were used by Kulikova et al. (2004) and Jin et al. (2014), which restrict their power to understand species’ evolutionary histories (Degnan and Rosenberg 2009). Third, Kulikova et al. (2004) explored the potential causes of non-monophyly by comparing genetic diversity at mitochondrial and nuclear loci; however, this method is arguably unreliable due to the stochasticity of the coalescent process (Hudson and Turelli 2003; Maddison and Knowles 2006; Peters et al. 2007). Fourth, although interspecific gene flow was inferred in Kulikova et al. (2004), the magnitude of gene flow has not yet been measured. Thus, studies that use broader geographical sampling, greater numbers of nuclear loci and advanced analytic approaches are critical to obtaining a comprehensive understanding of the evolutionary history and relationship of the 2 duck species.

In this study, we integrated data from the mitochondrial DNA (mtDNA) control region and 14 nuclear loci derived from a broad sampling of Chinese spot-billed ducks with published sequences from Old World Mallards (Kulikova et al. 2004; Peters et al. 2014). We aimed to investigate whether the previously reported non-monophyly between Chinese spot-billed ducks and Mallards is also evident when considering larger sample sizes and multiple nuclear loci. Then, we tested for incomplete lineage sorting versus introgression and explored the direction and extent of gene flow. Under an introgression scenario, we expect the coalescent-based isolation-with-migration (IM) model to identify the model with bidirectional/unidirectional gene flow as the best-fit model, lineage sharing to be more common in sympatric areas than in allopatric areas, and divergence time larger than the estimated complete lineage sorting time required. Under an incomplete lineage sorting scenario, we expect the IM model to identify the strict allopatric speciation model as the best-fit model, lineage sharing to be randomly distributed across the ranges of the 2 species, and divergence time smaller than the estimated complete lineage sorting time required.

Materials and Methods

Sample collection and laboratory methods

We collected 83 samples from 16 sites spanning the major geographical distribution range of Chinese spot-billed ducks (Figure 1, Supplementary Table S1). Samples contained feathers and toepads along with some muscle, blood, and egg membranes. Feathers and egg membranes were collected from bird nests during the breeding season (April–June) in 2015. Toepads were provided by the National Zoological Museum of China. Genomic DNA from muscle, blood, and egg membranes was extracted using the DNeasy Tissue Kit (Qiagen, Hilden, Germany) and following the manufacturer’s protocol, and genomic DNA from feathers and toepads was extracted following a modified protocol (see Irestedt et al. 2006). We amplified the mtDNA control region, one Z-linked nuclear locus Chromo-helicase-DNA binding protein gene 1 (CHD1Z) and 13 autosomal loci (see Supplementary Table S2; Kulikova et al. 2004; Peters et al. 2012). For each sample, we sequenced the mtDNA control region, Ornithine decarboxylase intron 6 (ODC6), Serum amyloid A (SAA). For samples with high DNA template quality (muscle, blood, feathers, and egg membranes), we subsampled 5 individuals from each sampling site and sequenced the remaining 12 nuclear loci. In total, 31 samples from 7 sampling sites covering the region from southwestern to northeastern China (5 samples from Aer Mountain (AE), Xianghai (XH), Dongting Lake (DT), Poyang Lake (PY), Shengjin Lake (SJ), and Huize (HZ), and one sample from Huayin (HY)) were sequenced for all loci. Primers and protocols from the studies of Kulikova et al. (2004) and Peters et al. (2012) were followed. Complete sequences were assembled using SEQMAN II (DNASTAR) and compared visually to the original chromatograms to avoid reading errors. Sequences of Chinese spot-billed ducks from Primorye and sequences of Old World Mallards were obtained from published datasets (Kulikova et al. 2004; Peters et al. 2014).

Genetic diversity and relationships among individuals

Sequences were aligned using CLUSTALW implemented in MEGA 6 (Tamura et al. 2013). We phased heterozygous nuclear sequences using the program PHASE (Stephens et al. 2001). Only individuals with resulting phase probabilities greater than 0.7 were used in the subsequent analyses (Carling et al. 2010). Fu and Li’s D test and Tajima’s D test implemented in DnaSP 5.0 (Librado and Rozas 2009) were used to examine the selective neutrality of the mtDNA control region and the nuclear loci. DnaSP was also used to calculate genetic diversity, including polymorphic sites, number of haplotypes, and haplotype and nucleotide diversity. Global nuclear DNA (nuDNA) diversity was calculated based on 7 loci (Annexin A11 (ANXA11), CHD1Z, Alpha-B crystallin (CRYAB), Fibrinogen beta chain (FGB), Growth hormone 1 (GH1), Preproghrelin (GHRL), Phosphenolpyruvate carboxykinase (PCK1)). Other loci were excluded from the diversity estimation because of limited numbers of sequences.

Relative divergence (FST) between Mallards and Chinese spot-billed ducks was calculated across markers in Arlequin 3.5 (Excoffier and Lischer 2010). Global nuDNA FST was also calculated. Next, relationships among samples were explored through median-joining networks reconstructed across loci in Network 5.0 (Bandelt et al. 1999). In addition, neighbor-net trees using uncorrected P-distances were reconstructed for mtDNA and concatenated, consensus nuclear data in SplitsTree 4.14.2 (Huson and Bryant 2006). TREEFINDER (Jobb et al. 2004) was used to concatenate consensus sequences. Lecithin-cholesterol acyltransferase (LCAT) and ODC6 were excluded from sequence concatenation because of the limited numbers of sequences. Heterozygous positions were coded with International Union of Pure and Applied Chemistry (IUPAC) ambiguity symbols and handled with average states in the SplitsTree analysis. Finally, we used a Bayesian approach implemented in STRUCTURE 2.3.4 (Falush et al. 2007) to assign individuals to clusters based on nuclear loci. We used xmfa2struct (available from http://www.xavierdidelot.xtreemhost.com/clonalframe.htm) to convert sequence data to STRUCTURE compatible files. ODC6 was deleted from the STRUCTURE analysis because only sequences from Primorye were available for Mallards. We used an admixture model without prior information regarding sampling locations. Analyses were run using correlated allele frequencies for 500,000 generations after a burn-in of 500,000 steps. The number of groups (K) was varied from 1 to 8, with 5 runs for each K value. The optimum K was determined by calculating ΔK in the program STRUCTURE HARVESTER (Earl 2012). CLUMPP (Jakobsson and Rosenberg 2007) was used to find the optimal alignments of the 5 replicate runs using the “FullSearch” algorithm. DISTRUCT 1.1 (Rosenberg 2004) was used to display the results.

IM analysis

To understand the demographic histories of Chinese spot-billed Duck and Mallard, we applied the coalescent-based IM model implemented in the program IMa (Hey and Nielsen 2004; Hey 2005). The IM model simultaneously estimates 6 parameters scaled to the mutation rate: ΘC (effective population size of Chinese spot-billed ducks), ΘM (effective population size of Mallards), ΘA (effective population size of the ancestral population at the time of divergence), t (time since divergence), m (migration rate from Mallards to Chinese spot-billed ducks), and m (migration rate from Chinese spot-billed ducks to Mallards). We used IMgc (Woerner et al. 2007) to filter for intralocus recombination using the default parameter setting.

To facilitate comparisons of gene flow between different marker systems, we conducted IM analyses based on the mtDNA and nuDNA datasets separately. All of the analyses were initially run in “M mode”. We applied the infinite site (IS) model of nucleotide substitution to the nuclear loci and the Hasegawa, Kishino, and Yano (HKY) model to the mtDNA control region as suggested by Hey (2007). The inheritance scalar was set to 1 for the autosomal loci, 0.75 for the Z-linked locus and 0.25 for the mtDNA control region to account for differences in effective population size. To improve mixing, a geometric heating scheme with 6 chains was used. We performed multiple runs, with an increasing number of steps and wide priors for each parameter. Based on these runs, we defined narrower priors that encompassed the full posterior distribution of each parameter. We assessed the convergence of the Markov chain Monte Carlo (MCMC) chain by monitoring trend plots and requiring all effective sample sizes (ESS) to be >150 (Peters et al. 2007). We also repeated the analysis 3 times with a different random number seed to confirm that the replicates converged on the same stationary distributions. For the mtDNA dataset, analyses were run for 30 million steps with a burn-in period of 3 million steps. For the nuDNA dataset, estimates of the migration rate contained several peaks, and the upper bounds did not return to zero regardless of prior adjustment and running step increase. Limited genetic variation might prevent the nuDNA gene flow parameters from being estimated accurately. After the M mode runs, we further ran IMa in “L mode” to compare likelihoods among 4 models with free or constrained migrations: model with bidirectional gene flow (ABCDD), model with unidirectional gene flow from Mallards to Chinese spot-billed ducks (ABCD0), model with unidirectional gene flow from Chinese spot-billed ducks to Mallards (ABC0D), and strict allopatric speciation model (ABC00).

Results

mtDNA nucleotide diversity of each population ranged from 0.004 to 0.017 in Chinese spot-billed ducks and from 0.004 to 0.010 in Mallards, and haplotype diversity ranged from 0.55 to 1 in Chinese spot-billed ducks and from 0.86 to 1 in Mallards (Supplementary Table S1). Nucleotide diversity across nuDNA loci ranged from 0.001 to 0.021 in Chinese spot-billed ducks and from 0.000 to 0.020 in Mallards, and haplotype diversity ranged from 0.16 to 0.98 in Chinese spot-billed ducks and from 0.04 to 0.98 in Mallards (Supplementary Table S2). The global nucleotide diversity and haplotype diversity of nuDNA were similar between the 2 species. Tajima’s D test revealed that none of the loci departed significantly from neutrality except the control region in Mallards and CRYAB in Chinese spot-billed ducks. The assumption of neutrality could not be rejected by Fu and Li’s D test for any loci except Nucleolin (NCL) in Chinese spot-billed ducks.

In the 545-bp mtDNA alignment, 74 sites were variable and 55 were parsimony informative. Polymorphic sites defined 107 haplotypes, of which only 5 haplotypes were shared between the 2 species. The mitochondrial haplotype network (Supplementary Figure S1) and neighbor-net tree (Figure 2A) recovered 2 divergent clades separated by 12 mutation steps, corresponding to Kulikova et al.’s (2004) group A and group SB. Group A was represented by 96 haplotypes, in which Mallards (61.1%, N = 118) were found at higher proportions than Chinese spot-billed ducks (38.9%, N = 75). Sixty-six haplotypes from group A were only observed in Mallards, and 25 were only observed in Chinese spot-billed ducks. Group SB was represented by 11 haplotypes, in which Chinese spot-billed ducks (87.9%, N = 29) were found at higher proportions than Mallards (12.1%, N = 4). Eight haplotypes from group SB were only observed in Chinese spot-billed ducks, and 3 were only observed in Mallards. In general, Mallards (96.7%) and Chinese spot-billed ducks (72.1%) tended to harbor group A haplotypes; however, the global FST value between the 2 species was 0.13, likely driven by the higher occurrence of the group SB in Chinese spot-billed ducks (27.9%). The FST reduced to 0.07 if the group SB was excluded. Mitochondrial haplotype distribution was not geographically structured (Figure 1).

Figure 2.

Neighbor-net trees of (A) mtDNA and (B) concatenated nuDNA loci. For the mtDNA tree, black circles represent haplotypes from Chinese spot-billed ducks, white circles represent haplotypes from Mallards, and gray circles represent haplotypes from both species. For the nuDNA tree, black circles represent sequences of Chinese spot-billed ducks, and white circles represent sequences of Mallards.

Whereas the 2 species showed low levels of divergence across the autosomal markers (average FST = 0.014; range −0.02 to 0.07; Supplementary Table S2), the Z-linked locus CHD1Z yielded a ten-fold higher FST estimate (FST = 0.15). With the exception of the Z-linked locus CHD1Z, nuDNA haplotypes provided no observable structuring between Mallards and Chinese spot-billed ducks across the remaining 11 loci (Supplementary Figure S2). No structure was further observed in the reconstructed neighbor-net tree (Figure 2B) or in STRUCTURE (Figure 3). Based on ΔK values, the optimal number of clusters was 4 (Supplementary Figure S3). However, the selecting K does not make clear assignments either in species or geographical distribution, and this result should be treated with caution because of no clear biological interpretation (Pritchard et al. 2010). At K = 2, 15 of the 23 Mallards and 7 of the 25 Chinese spot-billed ducks were assigned to cluster one with an average probability of 87.7% (±9.4% SD), and the remaining individuals were assigned to cluster two with an average probability of 91.0% (±10.5% SD). Cluster two was randomly distributed in the 2 species. Cluster one were randomly distributed in Mallards, but were restricted in northeastern (AE and XH) and southern China (HZ) in Chinese spot-billed ducks. At K = 3, 7 Mallards were further assigned to cluster three with an average probability of 96.6% (±1.4% SD). Higher values of K did not provide additional and interpretable resolution of the population structure.

Figure 3.

STRUCTURE assignment probabilities based on the nuDNA dataset (K = 2–4). Each vertical bar represents an individual’s assignment to a genotypic cluster with colors designating the different clusters.

The IMa analyses based on the mtDNA dataset produced posterior distributions of all parameters with clear peaks and bounds within the prior distributions (Figure 4). The posterior distributions of ΘC peaked at 79.22 (90% highest posterior density [HPD]: 49.62–118.54), ΘM peaked at 402.2 (263.8–659.8 HPD), m peaked at 0.01 (0.01–0.50 HPD), m peaked at 0.25 (0.06–0.54 HPD), and t peaked at 1.04 (0.80–1.32 HPD). These estimates were converted into biologically meaningful values using an assumed substitution rate of 4.8 × 10−8 per site per year and a generation time of 3 years (Peters et al. 2005; Peters 2006; see Hey 2007 for details on conversions). The effective population sizes of Chinese spot-billed ducks and Mallards were estimated as 252,000 (158,000–377,000 HPD) and 1,281,000 (840,000–2,101,000 HPD), respectively. The migration rate from Chinese spot-billed ducks to Mallards was ∼49.67 (7.52–177.17 HPD) migrants per generation, while that from Mallards to Chinese spot-billed ducks was negligible, ∼0.04 (0.02–29.81 HPD) migrants per generation. Time since divergence between the 2 duck species was estimated to be 40,000 (30,000–50,000 HPD) years ago. Each of the 4 models with free or constrained migration was rejected by the IMa analyses (P < 0.001; Table 1). However, the model with unidirectional gene flow from Chinese spot-billed ducks to Mallards (ABC0D) was superior to the other 3 models. As the estimates of nuDNA migration rate were unreliable, we did not include the results of the IMa analyses based on the nuDNA dataset.

Figure 4.

Demographic parameters estimated with IMa based on the mtDNA dataset. All parameter estimates are scaled to mutation rate. ΘC, effective population size of Chinese spot-billed ducks; ΘM, effective population size of Mallards; m, migration rate from Mallards to Chinese spot-billed ducks; m, migration rate from Chinese spot-billed ducks to Mallards; t, time since divergence.

Table 1.

Summary of likelihood ratio test statistics from the nested model analysis in IMa

Model	2LLR	df	P
ABCDD (mC = mM > 0)	39.791	1	<0.001
ABCD0 (mC > 0, mM = 0)	94.148	1	<0.001
ABC0D (mC = 0, mM > 0)	29.386	1	<0.001
ABC00 (mC = mM = 0)	123.716	2	<0.001

mC, migration rate from mallards to Chinese spot-billed ducks; mM, migration rate from Chinese spot-billed ducks to mallards; 2LLR, 2*log likelihood ratio; df, degrees of freedom; P-value, significance level; ABCDD, model with bidirectional gene flow; ABCD0, model with unidirectional gene flow from mallards to Chinese spot-billed ducks; ABC0D, model with unidirectional gene flow from Chinese spot-billed ducks to mallards; ABC00, strict allopatric speciation model.

Discussion

Genetic differentiation between Chinese spot-billed ducks and mallards

To date, we provide the most geographically comprehensive study on the phylogeography of Chinese spot-billed ducks and sympatric Mallards. Although a number of nuclear and mitochondrial markers were used to explore population structure, we report little differentiation (Supplementary Table S2) and an inability to successfully assign samples to their respective taxon (Figures 1–3; Supplementary Figure S1). We note that STRUCTURE did assign Chinese spot-billed ducks from northeastern (AE and XH) and southern China (HZ) to a single, unique cluster. Alleles from northeastern China (sympatric areas with Mallards) might be introgressed from Mallards. However, alleles from HZ were unlikely to come from Mallards because Chinese spot-billed ducks from HZ are sedentary and allopatrically distributed with Mallards (Carboneras and Kirwan 2018a; del Hoyo et al. 2018). Further research on the phylogeography of Indian Spot-billed Ducks (A. poecilorhyncha), a species found breeding sympatrically and hybridizing with Chinese spot-billed ducks in southern China might shed light on the formation of the phylogeographical structure (Leader 2006; Carboneras and Kirwan 2018b).

Subtle interspecific divergence was recovered in the mtDNA control region and Z-linked locus CHD1Z, with the major allele frequencies differing between the 2 species (Figure 2A, Supplementary Figures S1 and S2). Interspecific differentiations were ∼10 times larger at the mtDNA control region and CHD1Z than at the autosomal loci. Both mtDNA and Z-linked markers are often found to have higher divergence than autosomal markers (Sætre et al. 2003; Hogner et al. 2012; Toews and Brelsford 2012; Dhami et al. 2016). The effective population sizes of mtDNA and Z-linked markers are 1/4 and 3/4 that of autosomal markers, respectively. Based on genetic drift alone, the divergences of mtDNA and Z-linked markers should be 4-fold and ∼1.33-fold larger than autosomal markers, respectively (Zink and Barrowclough 2008; Ellegren 2009; Lavretsky et al. 2015a). Thus, we conclude that the 10-fold differences observed in the FST ratios between mtDNA or CHD1Z and autosomal markers cannot be explained by genetic drift alone.

Two of the most commonly cited causes of discordance between mtDNA and autosomal markers are selection and sex-biased dispersal (Toews and Brelsford 2012; Peters et al. 2014). Despite the long-held assumption that mtDNA is a putatively neutral marker, a number of studies have suggested that mtDNA is under natural selection (Galtier et al. 2009; Toews and Brelsford 2012). Local adaptation could reduce gene flow and elevate mtDNA divergence (Ribeiro et al. 2011; Pavlova et al. 2013). That the control region in Mallards significantly departed from neutrality suggests that selection might contribute to the higher differentiation in mtDNA than in autosomal loci, though the departure from neutrality might also relate to population expansions. Haldane’s rule, which posits that hybrids of the heterogametic sex (female in birds) have lower fitness than the homogametic sex (Haldane 1992), may also reduce mtDNA gene flow and result in a high level of mtDNA differentiation. We do not think sex-biased dispersal contributed to the high mtDNA differentiation because ducks normally exhibit male-biased dispersal (Rohwer and Anderson 1988; Gay et al. 2004; Liu et al. 2012), which would decrease instead of increase interspecific mtDNA divergence (Petit and Excoffier 2009).

In birds, which have the ZW sex-determinant system, Z-linked markers are often found to have elevated levels of interspecific divergence compared with those of autosomal markers, especially during the early stages of speciation (Sætre et al. 2003; Hogner et al. 2012; Dhami et al. 2016). This scenario is because the Z chromosome has been demonstrated to play an important role in building up prezygotic (i.e., male plumage traits and female mating preferences) and postzygotic (i.e., low hybrid fitness) isolation between species (Sætre et al. 2003; Sæther et al. 2007; Storchová et al. 2010). Similar to our study, elevated CHD1Z divergence has been found within 3 ducks (Mallard, Eurasian wigeon A. penelope and green-winged teal A. crecca; Peters et al. 2014), between 2 flycatcher species (Ficedula hypoleuca and F. albicollis; Backström et al. 2010), and between 2 sparrow species (Passer domesticus and P. hispaniolensis; Elgvin et al. 2011). In addition, elevated Z-linked divergence has been documented between Mallards and Mexican Ducks (A. [platyrhynchos] diazi), which is likely the product of divergent selection (Lavretsky et al. 2015a). We think that the Z chromosome might also be playing an important role in the divergence of Chinese spot-billed ducks and Mallards, and that divergent selection has resulted in the elevated Z-linked differentiation observed in our study.

Incomplete lineage sorting and introgression in the formation of non-monophyly

Introgression has been suggested to be an important factor contributing to non-monophyly among many duck species (McCracken et al. 2001; Mank et al. 2004; Peters et al. 2007; Kraus et al. 2012; Lavretsky et al. 2015b). In some cases, the role of incomplete lineage sorting might be more important (Avise et al. 1990; Lavretsky et al. 2014a). Kulikova et al. (2004) indicated that introgression might be the principle cause of non-monophyly between Chinese spot-billed ducks and Mallards by comparing genetic diversity at mitochondrial and nuclear loci. Our study, integrating multiple analyses methods, indicates that both incomplete lineage sorting and introgression might contribute to the non-monophyly.

Several lines of evidence support the incomplete lineage sorting scenario. First, the 2 species showed extensive nuDNA sharing, which has been commonly interpreted as evidence of recent divergence and the retention of ancestral polymorphisms (Peters et al. 2007). Second, the random distributions of mtDNA group A and nuDNA cluster two are consistent with the incomplete lineage sorting scenario. Under an introgression scenario, lineage sharing might be more common in sympatric areas than in allopatric areas because species have more opportunities for hybridization in sympatric areas (Kulikova et al. 2005; Peters et al. 2007; Wang et al. 2014). In contrast, incomplete lineage sorting does not necessarily make any predictions regarding the geographic distribution of lineages. Third, the interspecific divergence time is much more recent than the estimated complete lineage sorting time required. For mtDNA, the transition to monophyly for recently diverged populations usually requires 1 × Nfe generations where Nfe is the effective female population size (Moore 1995). Our IMa analyses indicated that the effective population sizes of Chinese spot-billed ducks and Mallards were ∼252,000 and 1,281,000, respectively. Therefore, at least 252,000 years might be required for segregating alleles to sort among lineages. However, the estimated divergence time, ∼40,000 (30,000–50,000 HPD) years ago was much a younger split.

Introgression might also contribute to the non-monophyly between the 2 duck species. The IMa analyses based on the mtDNA dataset indicated the existence of interspecific gene flow (Figure 4 and Table 1). The mtDNA group SB was restricted to sympatric areas in Mallards, and the nuDNA cluster one was mostly restricted to sympatric areas in Chinese spot-billed ducks, supporting the introgression hypotheses. Furthermore, the introgression scenario is concordant with the well-documented evidence of hybridization between the 2 species (Melville 1997; Kulikova et al. 2004).

Evolutionary histories of the 2 species

Omland (1997) argued that the various monochromatic species in the Mallard complex might have originated from a dichromatic Mallard-like common ancestor. Similarly, Avise et al. (1990) concluded that the American Black Duck (A. rubripes) is a recent evolutionary derivative of a more broadly distributed Mallard-Black ancestor, and the paraphyly of Mallards and American Black Ducks was the result of incomplete lineage sorting. We hypothesize that the Chinese spot-billed duck might also arise from a dichromatic Mallard-like ancestor, and the divergence was too recent for ancestral polymorphisms to have sorted into reciprocally monophyletic lineages. During the diversification of Chinese spot-billed ducks and Mallards, introgression further caused gene mixture between the 2 species. The mtDNA gene flow was unidirectional from Chinese spot-billed ducks to Mallards (Figure 4). In agreement with the unidirectional gene flow, group SB haplotypes were widely distributed in Chinese spot-billed ducks but were restricted to sympatric areas in Mallards. Their high frequency in Chinese spot-billed ducks suggests that these haplotypes might have originated from Chinese spot-billed ducks and been introduced into Mallards via hybridization (Kulikova et al. 2005, 2012).

The asymmetric mtDNA introgression might result from reciprocal hybridization, but the offspring resulting from crosses of Mallard females and Chinese spot-billed duck males are inviable or infertile (asymmetric postzygotic isolation). Asymmetric postzygotic isolation is very common (Haldane 1922; Turelli and Moyle 2006). However, data on the fitness of hybrids between Chinese spot-billed ducks and Mallards are lacking; thus, we cannot assess the probability of this explanation. In addition, the asymmetric mtDNA introgression might result from unidirectional hybridization, with offspring formed mostly from Chinese spot-billed duck females and Mallard males (asymmetric prezygotic isolation). In general, females prefer to mate with conspecific males (Brodsky and Weatherhead 1984; Baker 1996) but sometimes mate with heterospecifics when the latter dominate in number (Avise and Saunders 1984; Wirtz 1999; McCracken and Wilson 2011), possess brighter plumage (Norris 1990; Stein and Uy 2006; Wang et al. 2014), or behave more aggressively (McDonald et al. 2001; Krosby and Rohwer 2009). Mallard males are more colorful in plumage than Chinese spot-billed duck males. Moreover, Mallard males are typically dominant over males of other duck species and are commonly successful when competing for females (Brodsky and Weatherhead 1984; Brodsky et al. 1988; but see Hoysak and Ankney 1996). Therefore, hybridization between Mallard and closely related species primarily involves Mallard males and females of other species (Brodsky and Weatherhead 1984; Brodsky et al. 1988; but see Fowler et al. 2009). Similarly, hybridization between Chinese spot-billed ducks and Mallards might also primarily involve Mallard males and Chinese spot-billed duck females.

In conclusion, applying widespread sampling, multiple nuclear loci, and coalescent-based analytic approaches, the present study confirms genetic sharing between 2 parapatric duck species in East Asia. Furthermore, both incomplete lineage sorting and introgression might contribute to the observed non-monophyly, implying complicated evolutionary history. Given that the 2 duck species are morphologically differentiated, but weakly differentiated at autosomal loci, we hypothesize that a few genes, probably from the Z chromosome, might play a major role in the divergence and maintenance of species integrity in the face of gene flow (Wu 2001; Feder et al. 2012). Genome-scale data are thus necessary to uncover genomic architecture and the associated underlying mechanisms and may provide further insights into the evolutionary histories of the 2 species and other waterfowls.

Click here for additional data file.