Skip to main content
Download PDF

Past hybridization between two East Asian long-tailed tits (Aegithalos bonvaloti and A. fuliginosus)

Frontiers in Zoology volume 11, Article number: 40 (2014) Cite this article

Abstract

Introduction

Incomplete lineage sorting and hybridization are two major nonexclusive causes of haplotype sharing between species. Distinguishing between these two processes is notoriously difficult as they can generate similar genetic signatures. Previous studies revealed that the mitochondrial DNA (mtDNA) differentiation between two East Asian long-tailed tits (Aegithalos bonvaloti and A. fuliginosus) was extremely low, even lower than intraspecific differentiation in some other long-tailed tits. Using a combination of multilocus and coalescent analyses, we explored the causes of the anomalous lack of mtDNA differentiation between the two species.

Results

The mtDNA divergence between the two species was shallow, while the nuclear DNA (nuDNA) divergence was considerably deeper. The IMa analyses based on the mtDNA dataset suggested relatively high gene flow from A. fuliginosus to A. bonvaloti, while negligible gene flow in the opposite direction. In contrast to mtDNA, the migration rates at autosomal and Z-linked nuDNA loci were negligible or much lower. The NEWHYBRIDS analysis assigned all individuals except one to pure parental species with high posterior probability. The Bayesian skyline plot showed that both species underwent population expansions during the Last Glacial Maximum (LGM), and the ecological niche modelling suggested that their ranges overlapped more during the LGM than at present.

Conclusions

We suggest that historical hybridization, in combination with selective sweep and/or genetic drift might be the main causes of the extremely low mtDNA differentiation between the two species. The hybridization probably occurred mainly between A. fuliginosus females and A. bonvaloti males. The LGM distribution expansion might have facilitated hybridization, while the post-LGM distribution contraction could have facilitated some mtDNA sorting. Ongoing hybridization between the two species might be very limited, but further studies with more samples from the contact zone are needed to test this conclusion.

Introduction

It is widely recognized that gene trees can differ from the species phylogeny, and gene genealogies usually proceed through various stages of non-monophyly during the early stages of their divergence history[1–6]. In a survey of published mitochondrial analyses, Funk & Omland[4] identified species-level non-monophyly in 23% of the more than 2300 assayed species. Hybridization and incomplete lineage sorting (ILS) are considered to be two of the major, nonexclusive causes of species-level non-monophyly[4, 7, 8]. Hybridization between species may result in the incorporation of alleles from one species into the gene pool of another[9]. When the first generation hybrids backcross with at least one of the parental species, introgression occurs. Alternatively, ILS may also result in species-level non-monophyly if species divergence was too recent for ancestral polymorphisms to have sorted to reciprocal monophyly. Distinguishing between hybridization and ILS is notoriously difficult, because they leave similar genetic signatures[10–12].

A commonly used approach to distinguish between hybridization and ILS is comparing genetic diversity between mitochondrial DNA (mtDNA) and nuclear DNA (nuDNA). The presence of fixed differences between species in nuDNA despite mtDNA haplotype sharing, or deeper divergences in nuDNA than in mtDNA, are often interpreted as evidence of mtDNA migration[13–17]. However, the stochasticity of the coalescent process (i.e. randomness of mutation and genetic drift) and the potential for locus-specific effects such as selection or recombination, both of which might result in fixed base pair differences between species at some loci while shared haplotypes at other loci, prevent strong conclusions[18]. Massive hybridization over a long period could also homogenize nuDNA, making conclusions more difficult to draw[19]. Most importantly, this approach cannot exclude the possibility of a combined effect of both hybridization and ILS. Therefore, quantitative evaluation of the likelihood of alternative hypotheses is needed to draw robust conclusions[15, 20].

The development of the isolation-with-migration (IM) coalescent model, which takes into account inheritance models of markers and the stochasticity of the coalescent process, allows us to rigorously test explicit hypotheses[21–23]. The IM model is usually used to test two hypotheses related to divergence processes. The first hypothesis is a strict isolation model whereby speciation resulted from allopatric divergence without subsequent gene flow. The second hypothesis is an isolation with migration model that allows speciation despite the presence of gene flow. Non-monophyly can be attributed to ILS when the IM model supports the strict isolation model. In contrast, if the IM model supports the isolation with migration model, introgression may be necessary to explain non-monophyly. Additionally, the IM model could also be used to reveal introgression rates and gene flow directions at different marker systems[24–26], and identify candidate gene for reproductive isolation[27, 28].

Aegithalos bonvaloti (Black-browed Tit) and A. fuliginosus (Sooty Tit) are two East Asian endemic long-tailed tits. The two species are closely related to each other[29, 30]. A. bonvaloti is distributed from southwest China to extreme northeast Myanmar, and A. fuliginosus is restricted to the mountainous areas of central China (Figure 1;[31, 32]). The two species are parapatrically distributed, with distributions partially overlapping to the west of the Sichuan Basin. A. bonvaloti is distributed from about 1500 – 4400 m, and A. fuliginosus from about 1000 – 2600 m[33]. Two subspecies are recognized in A. bonvaloti: bonvaloti in the southern part of the species’ range and a doubtful subspecies obscuratus (see Discussion for details) in the northern part[30]. In contrast, A. fuliginosus is monotypic. The plumage differences between the two species are pronounced (cf. Figure 1). Intriguingly, the mtDNA divergence between them is extremely low (cytochrome-b uncorrected p-distance: 0.002), even lower than intraspecific differentiation within some other long-tailed tit species[29, 30]. Päckert et al.[30] suggested that A. bonvaloti and A. fuliginosus may “represent just distinct morphotypes at an intraspecific differentiation level”. However, both Dai et al.[29] and Päckert et al.[30] included only one sample of A. bonvaloti and fewer than six samples of A. fuliginosus, so whether their results represent the true level of differentiation is unknown.

Figure 1
figure 1

Map of the study area and sampling locations. Solid line represents the distribution of Aegithalos bonvaloti and A. fuliginosus. Dotted line represents the approximate contact zone. Circles and triangles show the sampling locations of A. bonvaloti and A. fuliginosus, respectively. The colours of circles and triangles correspond to mitochondrial DNA (mtDNA) origin, with white for A. fuliginosus haplotypes and black for A. bonvaloti haplotypes. Elevation is represented by shading: lower elevations are shown lighter and higher elevations are shown darker. Insert: The dotted box represents the study region.

Full size image

In the present study, we explored mtDNA diversification between A. bonvaloti and A. fuliginosus using multiple samples throughout the range of both species. We first investigated whether the previously reported extremely low mtDNA divergence between them is also evident under larger sample sizes. As this was confirmed, we compared the genetic variation between mtDNA and nuDNA datasets, and utilized the IM model to test whether this was caused by ILS or hybridization. Finding support for historical hybridization, we tried to explore the interspecific gene flow pattern and related explanations. Finally, as it is uncertain whether the individuals representing intermediate phenotypes from the contact zone are of hybrid origin or represent a distinct subspecies A. b. obsuratus, we further investigate if contemporary hybridization occurs between the two species.

Results

Genetic polymorphism and species phylogeny

A total of 26 and 31 individuals were sampled for A. bonvaloti and A. fuliginosus, respectively (Figure 1 and Table 1). Two mtDNA fragments, three autosomal nuDNA loci and three Z-linked nuDNA loci were utilized. The mtDNA fragments were cytb (cytochrome b) and ND2 (NADH dehydrogenase subunit 2). The three autosomal nuDNA loci were MB (myoglobin, intron 2), TGFB2 (transforming growth factor beta-2, intron 5) and RAG1 (recombination activating protein 1). The three Z-linked nuDNA loci were BRM15 (brahma protein, intron 15), CHDa (chromo-helicase-DNA-binding protein, intron a) and VLDLR9 (very low density lipoprotein receptor, intron 9).

Table 1 Mitochondrial DNA nucleotide polymorphism at each sampling site
Full size table

The absence of ambiguous sites and stop codons suggested that these mtDNA sequences were not “numts”[34]. The combined mtDNA alignment contained 21 polymorphic sites, which defined 14 haplotypes. The mtDNA haplotype diversity of each population ranged from 0.40 to 1.00 in A. bonvaloti and from 0.68 to 1.00 in A. fuliginosus, and the nucleotide diversity of each population ranged from 0.00041 to 0.00157 in A. bonvaloti and from 0 to 0.0021 in A. fuliginosus (Table 1). For the nuDNA loci, the number of haplotypes, haplotype diversities and nucleotide diversities were similar between the two species (Table 2). The results of multilocus HKA test demonstrated no significant departure from neutrality for all loci [sum of deviations = 18.40, degrees of freedom (df) = 12, P = 0.10]. Tajima’s D test[35] revealed that most loci did not depart significantly from neutrality, except for ND2 and the entire mtDNA alignment in A. fuliginosus. The assumption of neutrality for all loci could not be rejected by Fu and Li’s D test[36].

Table 2 Summary statistics for each gene region and values of neutrality tests
Full size table

The two species were not reciprocally monophyletic in the mtDNA gene tree (Figure 2). The most abundant mtDNA haplotype was shared by 19 individuals of A. fuliginosus and 5 individuals of A. bonvaloti. No other mitochondrial haplotype was shared between the two species. Nuclear haplotypes were shared between the two species at BRM15, VLDLR9 and RAG1 (Figure 3). However, fixed base pair differences were found at MB, TGFB2 and CHDa.

Figure 2
figure 2

Maximum clade credibility Bayesian tree and haplotype network of mitochondrial DNA. Branch support values higher than 0.80 are shown above nodes. Each circle represents a haplotype, and the size of the circle is proportional to that haplotype’s frequency. Dots represent unsampled haplotypes. Haplotypes are coloured according to taxon: black for Aegithalos bonvaloti and white for A. fuliginosus. H2 is the most widely distributed haplotype shared by 19 individuals of A. fuliginosus and 5 individuals of A. bonvaloti.

Full size image
Figure 3
figure 3

Haplotype networks constructed with nuclear loci. The three loci in the upper row (CHDa, VLDLR9 and BRM15) are Z-linked, whereas the three in the lower row (TGFB2, RAG1 and MB) are autosomal. Each circle represents a haplotype, and the size of the circle is proportional to that haplotype’s frequency. Dots represent unsampled haplotypes. Haplotypes are coloured according to taxon: black for Aegithalos bonvaloti and white for A. fuliginosus.

Full size image

In order to evaluate whether the mtDNA divergence is reduced relative to nuDNA, we calculated the uncorrected p-distance between the two species at each locus, and the mtDNA to nuDNA divergence ratios between the two species and other Aegithalos species pairs. The mean uncorrected p-distance between A. bonvaloti and A. fuliginosus was 0.002 for mtDNA, whereas it was on average twice as high for the nuDNA loci (mean 0.004; range 0.001 to 0.008; Table 2). The mean ratio of mtDNA to nuDNA divergence of A. bonvaloti and A. fuliginosus (0.4) was considerably lower than in other Aegithalos species pairs (mean 7.6; range 4.4 to 10.7; Additional file1).

Coalescent estimates of gene flow

We applied the coalescent-based IM model implemented in the program IMa[22] to distinguish between ILS and hybridization. This method co-estimates several parameters, including effective population size, divergence time and migration rate. In order to reduce the risk that the IMa would confuse genetic homoplasy with hybridization, we conducted two kinds of analyses, one with all individuals, and the other with only allopatric individuals. If the analysis based on all individuals reveals interspecific gene flow, while the analysis based on only allopatric individuals reveals near-zero gene flow, this suggests that the interspecific gene flow is restricted to the hybrid zone, and that hybridization occurs between the species. However, if both analyses show interspecific gene flow, this suggests that ILS might be needed to explain the genetic pattern. To facilitate comparisons of introgression estimates between different marker systems, we conducted analyses based on three separate datasets: one with the three Z-linked loci, one with the three autosomal loci, and one with the mtDNA.

The IMa analyses produced posterior distributions of migration rate and effective population size parameters with clear peaks and bounds within the prior distribution (Figure 4 and Additional file2), which suggests that these results were reliable. The posterior distributions of the divergence time parameters also had clear peaks, however the upper bounds did not return to zero within the prior distribution. The low levels of variation in the analyzed loci may prevent the divergence time parameters from being estimated accurately.

Figure 4
figure 4

The marginal posterior probability density distributions of migration rates ( m ). All parameter estimates are scaled to the mutation rate. The black line represents migration rate from Aegithalos fuliginosus to A. bonvaloti. The grey line represents migration rate from A. bonvaloti to A. fuliginosus. “all” and “allopatric” mean IM analyses based on all individuals and allopatric individuals, respectively.

Full size image

For the mtDNA dataset, none of the four models could be rejected (P> 0.05, Additional file3). For the autosomal and Z-linked DNA datasets, the strict allopatric speciation models were rejected. For the mtDNA dataset including all individuals, the estimated migration rate from A. fuliginosus to A. bonvaloti was about 1.54 migrants per generation, while the migration rate from A. bonvaloti to A. fuliginosus was negligible (Table 3). For the mtDNA dataset based on allopatric individuals, the estimated gene flow in both directions was negligible. For the autosomal DNA dataset, the migration rates in both directions were near zero. For the Z-linked DNA dataset including all individuals, the migration rate from A. bonvaloti to A. fuliginosus was about 0.18 migrants per generation, while the migration rate from A. fuliginosus to A. bonvaloti was negligible. For the Z-linked DNA dataset based on allopatric individuals, the estimated migration rate from A. bonvaloti to A. fuliginosus was 0.15 migrants per generation, i.e. similar to that based on all individuals.

Table 3 Demographic parameters estimated with IMa
Full size table

Historical demographic and spatial dynamics

To gain insight into the relationship between hybridization and past population dynamics, we explored the historical demographic and spatial dynamics of the two species. Mismatch distribution and Bayesian skyline plot (BSP) analyses were conducted to explore the past population dynamics. The mismatch distributions for both species were unimodal, which suggest that the two species underwent population expansions. The results of the BSP analysis suggested population expansions in both species during the Last Glacial Maximum (LGM; Figure 5).

Figure 5
figure 5

Bayesian skyline plots representing historical demographic trends of Aegithalos bonvaloti and A. fuliginosus . The x-axis is in unit of million years ago. The y-axis shows the scaled effective population size. The median estimate is enclosed within the 95% HPD limits. Black and grey represent A. bonvaloti and A. fuliginosus, respectively. Dotted lines show the Last Glacial Maximum (LGM).

Full size image

Ecological niche modelling (ENM) was used to predict the potential distributions of both species at present and at the LGM. The high AUC (area under the receiver operating characteristic curves; 0.991 for A. fuliginosus and 0.989 for A. bonvaloti) and low binomial probabilities (P≪ 0.0001) indicate that the ENMs for both species perform better than random predictions. The predicted present distributions are largely congruent with the species’ known distributions (Figure 6). During the LGM, the predicted distributions for both species were extended to the low elevational Sichuan Basin, resulting in potentially larger contact area between the two species. After the LGM, both species moved to high elevational areas, thus contracting the contact area.

Figure 6
figure 6

Ecological niche models for Aegithalos bonvaloti and A. fuliginosus . The LGM represents the Last Glacial Maximum. Circles and triangles represent the localities used to build ecological niche models for A. bonvaloti and A. fuliginosus, respectively. Blue represents the potential distribution for A. bonvaloti, red for A. fuliginosus, and dark blue for the contact zone between the two species. The dotted area shows the Sichuan Basin.

Full size image

Contemporary gene exchange

We used NEWHYBRIDS[37] to explore whether ongoing hybridization occurs between the two species. The NEWHYBRIDS analysis was based on the nuDNA dataset. The analysis assigned almost all individuals to pure parental species with high posterior probability (P> 0.99; Figure 7). Only one individual, an A. bonvaloti which had a foreign mtDNA haplotype, was assigned to pure A. bonvaloti with relatively low posterior probability (P = 0.78).

Figure 7
figure 7

Detection of ongoing hybridizaiton using NEWHYBRIDS. Each histogram represents an individual. Field identification of the individuals is labeled at the bottom of each histogram. The posterior probability for each individual to be assigned as pure Aegithalos bonvaloti, pure A. fuliginosus, F1 hybrid, F2 hybrid, backcross with A. bonvaloti (B × B), and backcross with A. fuliginosus (B × F) are labeled in different colors. Asterisks indicate A. bonvaloti individuals with foreign mtDNA haplotype.

Full size image

Discussion

Extremely low mtDNA divergence

Our results confirm the previous reports[29, 30] of extremely low mtDNA differentiation between A. bonvaloti and A. fuliginosus, and show that this low divergence prevails across a large geographical area. The two species were not reciprocally monophyletic in the mtDNA gene tree, and there were no fixed mtDNA differences between them. Five individuals of A. bonvaloti from sampling sites SM and BX shared a haplotype with individuals of A. fuliginosus from a wide geographical range. The uncorrected p-distance between the two species was lower than intraspecific divergence in most other long-tailed tit species[29, 30]. Moreover, the ratio of mtDNA to nuDNA divergence between the two species was much lower than in the other Aegithalos species pairs.

Several lines of evidence support that historical hybridization rather than ILS is the main cause of the exceptionally low mtDNA divergence between the two species. The fixed differences between the two species at MB, TGFB2 and CHDa suggest that they have experienced a prolonged period of independent evolution. This is in agreement with the significant morphological differences between them. Although the shallow divergence and lack of reciprocal monophyly in mtDNA is compatible with recent lineage splitting, that is unlikely given the considerably deeper divergences in nuDNA. The effective population size of mtDNA is one quarter of that of nuDNA, and the mutation rate of mtDNA is generally an order of magnitude greater than that of nuDNA, thus, the mtDNA should be considerably more differentiated than the nuDNA markers[38–40]. However, this is opposite to the observed pattern in the two species. The results of the IMa analyses provide additional quantitative support for hybridization. Although the IMa analysis based on the mtDNA dataset rejected neither the strict allopatric speciation model nor the models with gene flow, migration was detected between the two species. The IMa analysis based on all individuals showed high levels of mtDNA gene flow from A. fuliginosus to A. bonvaloti, while the analysis based on allopatric individuals showed near-zero gene flow. These results suggest that the estimated gene flow is restricted to the contact zone, and that hybridization has occurred between the two species. As further support, ILS would not be expected to result in haplotype sharing exclusively in a geographical overlap zone between the two species[4, 15, 41].

Our IMa analyses suggested that the mtDNA gene flow was asymmetric, with high gene flow from A. fuliginosus to A. bonvaloti, and negligible gene flow in the reverse direction. In agreement with this, we found an A. fuliginosus mtDNA haplotype at the northernmost sampling sites within the range of A. bonvaloti (BX and SM), whereas only native mtDNA haplotypes among all the samples of A. fuliginosus. These results suggest that hybridization occurred mainly between A. fuliginosus females and A. bonvaloti males.

In general, females would prefer to mate with conspecific males rather than with heterospecific ones when given a choice, but may sometimes mate with a heterospecific male when heterospecifics dominate in numbers[42, 43], possess brighter coloration[44], display more attractive courtship[45], or behave more aggressively[46, 47]. A female A. fuliginosus in an area with predominantly A. bonvaloti (i.e. the southern part of the distribution of A. fuliginosus) might mate with a male A. bonvaloti as A. bonvaloti is more brightly coloured than A. fuliginosus[31, 32], whereas a female A. bonvaloti in an area with mainly A. fuliginosus (i.e. the northern part of the distribution of A. bonvaloti) might be less likely to choose a male A. fuliginosus because of his duller coloration. Female-biased dispersal could also have contributed to the observed pattern, if females of A. fuliginosus have been more dispersive than males and more frequently dispersed southward to areas where conspecific males are uncommon. Northward dispersal of male A. bonvaloti could result in a similar pattern, but as the general dispersal pattern in birds is female-biased[48, 49], female-biased southward dispersal of A. fuliginosus seems more likely.

In contrast to mtDNA, the nuDNA gene flow was limited. Strong mtDNA introgression despite limited nuDNA migration has been found in several hybrid zones[8, 50–53]. This is consistent with the theoretical expectation that mtDNA is more likely than nuDNA to be transmitted across species boundaries[4, 54]. As the effective population size of mtDNA is one quarter of that of nuDNA, foreign mtDNA haplotypes stand a greater chance of reaching fixation than foreign nuDNA haplotypes[4]. Moreover, mtDNA is less likely to be constrained by linkage to loci under selection than nuDNA, and is thus more susceptible to move across species boundaries[55, 56].

However, we note that introgression cannot on its own explain the low mtDNA divergence between A. bonvaloti and A. fuliginosus. In view of the pronounced nuDNA divergences, the original mtDNA haplotypes should have been considerably more divergent than observed, resulting in sharing of a mixture of divergent haplotypes. As that is not the case, some process must have homogenized the mtDNA haplotypes. We suggest that a past selective sweep might have been involved. Given its role in oxidative phosphorylation, mtDNA is particularly likely to be the subject of selective sweeps, especially in organisms living in different thermal regimes[8, 57–59]. The two long-tailed tits indeed occur at partly different elevations with different environments. In addition, Tajima’s D test revealed that the ND2 gene and the entire mtDNA alignment in A. fuliginosus may be under natural selection (P< 0.05). In accordance with our results, natural selection on the ND2 gene has been found in human populations[60], and has been suggested to be the source of shallow genetic structure in Willow Tit (Parus montanus;[61]). A very similar pattern to that found in the two long-tailed tits, which was also suggested to be the result of a selective sweep, has been observed in two widespread Palearctic buntings, Yellowhammer Emberiza citrinella and Pine Bunting E. Leucocephalos[13, 17].

Genetic drift could also have increased the frequency of introgressed haplotypes, even though the introgressed haplotypes are not under positive selection[51, 62, 63]. Our BSP analyses revealed that the population sizes of both species were small before the LGM, which might have facilitated genetic drift.

In addition, range expansion during the LGM might have facilitated hybridization. Montane species generally undergo movements along elevational gradients to track the environmental conditions they have adapted over time[64, 65], descending to lower elevational areas during glacials and becoming vicariantly isolated on different mountains during interglacials. Glacial-induced range expansion has been suggested to create great opportunities for the contact and hybridization between populations associated with different mountain systems[66, 67]. The ENM analysis showed that both species underwent range expansions to the low altitude Sichuan Basin during the LGM, resulting in more extensively overlapping distribution. The larger overlapping distribution might facilitate hybridization between the two species.

Although the mtDNA divergence between A. bonvaloti and A. fuliginosus is extremely low, if the most abundant haplotype H2 is excluded, there are fixed differences between them. This suggests that the two species have undergone more recent mtDNA sorting. Our ENM indicated that both species’ ranges contracted after the LGM, thereby reducing the horizontal contact area. The partly different elevational ranges of the two species further limited the vertical contact area. The relatively narrow overlapping distributions might have allowed only a narrow conduit for mtDNA to pass between the two species, thus facilitating sorting of mtDNA haplotypes.

Ongoing gene exchange

Hybridization between A. bonvaloti and A. fuliginosus has been postulated based on a series of specimens representing an intermediate phenotype[30, 68]. However, A. bonvaloti has a doubtful subspecies, A. b. obscuratus, which is distributed at the contact zone between A. bonvaloti and A. fuliginosus. A. b. obscuratus is slightly duller and browner than A. b. bonvaloti, and has been regarded as approaching A. fuliginosus, or to be intermediate between A. b. bonvaloti and A. fuliginosus[30, 31]. It is uncertain whether the individuals representing intermediate phenotypes are hybrids or belong to A. b. obscuratus. Our NEWHYBRIDS analysis assigned all individuals except one to pure parental species with high posterior probability, suggesting that the ongoing gene exchange between the two species might be limited. However, we note that we included few samples from the contact zone, so more samples are needed to explore if hybridization still takes place, and to re-evaluate the status of A. b. obscuratus.

Conclusions

Based on larger sample sizes than in earlier studies, we confirm that the mtDNA divergence between A. fuliginosus and A. bonvaloti is extremely low. We suggest that past hybridization in combination with a selective sweep and/or genetic drift are the main causes of this pattern. The interspecific mtDNA gene flow was species-biased, with relatively high gene flow from A. fuliginosus to A. bonvaloti, while only negligible gene flow in the reverse direction. Non-assortative mating and/or female-biased southward dispersal of A. fuliginosus might have contributed to the unidirectional mtDNA gene flow. In contrast to mtDNA, the interspecific nuDNA gene flow was negligible or much lower. Range expansion during the LGM might have facilitated hybridization, while post-LGM distribution contraction could have led to some sorting of mtDNA haplotypes. Our results suggest that the ongoing interspecific gene exchange might be very limited. However, more samples from the contact zone are needed to draw robust conclusions.

Methods

Sampling and DNA sequencing

Sample sizes for A. bonvaloti and A. fuliginosus were 26 and 31 individuals, respectively. A. caudatus (Long-tailed Tit), which is closely related to A. bonvaloti and A. fuliginosus[29, 30], was used as outgroup. Samples were stored in 100% ethanol in the field and later at -80°C in the laboratory. Total genomic DNA was extracted from blood or muscle samples using the QIAamp DNA Mini Kit (QIAGEN) following the manufacturer’s protocol.

We amplified the two mtDNA fragments and six nuDNA loci via polymerase chain reaction (PCR). Primer pairs for cytb and ND2 were H16065/L14990 and H6313/L5219, respectively[69]. PCR conditions were as follows: denaturation at 94°C for 2 min, followed by 40 cycles at 93°C for 1 min, annealing temperature (49°C for cytb and 53°C for ND2) for 45 s, and 72°C for 1 min, and a final 8 min at 72°C. The primer pairs used to amplify the six nuDNA loci were Myo3F/Myo2, TGFB2.5F/TGFB2.6R[70], R50/R51[71], BRM15F/BRM15R[72], F2550/R2718[73], and VLDLR9F/VLDLR9R[72]. PCR conditions were as follows: an initial denaturation at 94°C for 2 min, 35 cycles of 94°C for 30 s, annealing temperature (51°C for BRM15, 61°C for RAG1 and 56°C for MB, TGFB2, CHDa and VLDLR9) for 30 s and 72°C for 30 s, and a final extension at 72°C for 10 min. Sequencing was carried out using the BigDye Terminator method (Applied Biosystems) on an ABI 3730xl analyzer. Complete sequences were assembled using Seqman II (DNASTAR) and compared visually to the original chromatograms to avoid reading errors. Sequences were aligned using CLUSTALW implemented in MEGA 5[74]. To avoid amplifying mtDNA homologues from the nuclear genome (“numts”), we checked whether the sequences had ambiguous sites or stop codons. We phased heterozygous nuclear sequences using the program PHASE[75]. Only individuals with resulting phase probabilities greater than 0.7 were used in subsequent analyses[76, 77]. All sequences have been deposited in GenBank (accession numbers KJ790262-KJ790672).

Genetic diversity and phylogenetic analyses

We concatenated sequences of the two mtDNA fragments to a single sequence of 1908 base pairs. Fu and Li’ D test and Tajima’s D test implemented in DnaSP 5.0[78] were used to assess whether nucleotide polymorphisms deviated from expectations under neutral theory for each locus in each species. Additionally, the multilocus HKA test was performed to test the neutrality of all loci using the HKA program[79]. DnaSP was used to calculate polymorphic sites, number of haplotypes, haplotype diversities and nucleotide diversities.

We used TCS 1.21[80] to produce a maximum-parsimony network for each locus. The Bayesian algorithm implemented in BEAST 1.7.0[81] was used to reconstruct phylogenetic relationships among mtDNA haplotypes. We selected the best-fit model of nucleotide substitution using the program MrModeltest 2.3[82]. According to the Akaike information criterion (AIC), the best fit model for mtDNA was the HKY model[83]. The uncorrelated lognormal relaxed molecular clock and Yule process tree prior were chosen. Analysis was run for 20 million steps and sampled every 2,000 generations. We assessed convergence of the Markov chain Monte Carlo (MCMC) chain in the program Tracer 1.5[84]. The output trees file was summarized into a maximum clade credibility tree using TreeAnnotator of the BEAST package. FigTree 1.3.1[85] was used to visualize the result.

We used MEGA to calculate uncorrected p-distance between the two species at each locus. We also calculated the ratio of mtDNA to nuDNA divergence between the two species and other Aegithalos species pairs. Three other Aegithalos species (A. niveogularis, A. concinnus and A. caudatus) were used. We used two mtDNA fragments (cytb and ND2) and three nuDNA loci [FIB (beta-fibrinogen, intron 7), TGFB2 and ODC (ornithine decarboxylase, intron 6–7)]. These species and loci were selected due to the accessibility in GenBank.

Coalescent simulations

In order to distinguish between ILS and hybridization, we applied the coalescent-based IM model implemented in the program IMa[22]. The IM model involves several simplifying assumptions, such as no recombination within loci, mutation following the model applied, no population structure within species, and no genetic contribution from unsampled species[22, 86]. A recent study has demonstrated that parameter estimates of the IM model are generally quite robust to small to moderate violations of the IM model assumptions[87]. We detected the minimum number of recombination events of each nuDNA locus using the four-gamete test in DnaSP. When recombination was detected, we kept the longest non-recombining block of sequences for subsequent analyses.

We conducted six IMa analyses. All analyses were first run in ‘M mode’. The HKY mutation model was used. The inheritance scalars were set as 1 for the three autosomal loci, 0.75 for the three Z-linked loci and 0.25 for the mtDNA to account for differences in effective population sizes. To improve mixing, we used a geometric heating scheme with six chains. We first performed multiple runs, with an increasing number of steps and using wide priors to ensure that the complete posterior distribution could be obtained. We finally performed three independent runs of 30 million steps with a burn-in period of three million steps. Convergence was assessed by monitoring trend plots, requiring all effective sample size (ESS) to be greater than 50, and ensuring the similarity of posterior distributions from independent runs. After M mode runs, we further ran IMa in ‘L mode’ to compare models with free or constrained migration.

Demographic reconstruction

We conducted mismatch distribution analysis using Arlequin 3.5[88]. Only the mtDNA alignment was used in this analysis. Individuals representing foreign mtDNA haplotype were excluded from this analysis. We also used the BSP method implemented in BEAST to estimate the posterior distributions of effective population size through time. Chains were run for 40 million generations, with the first 10% discarded as burn-in. The HKY model was used according to the result of MrModeltest. We used the uncorrelated lognormal relaxed molecular clock to account for rate variation among lineages. Demographic history through time was reconstructed with Tracer.

Ecological niche modelling

We used MAXENT[89] to predict the potential distributions of both species at present and at the LGM. We downloaded 19 bioclimatic variables at 2.5 arc-minute resolution from WorldClim ([90]; http://www.worldclim.org/). Species occurrence data included sampling sites, museum records (provided by the National Zoological Museum of China) and birdwatching locations (downloaded from http://birdtalker.net/report/). We removed sites separated from each other by less than 0.1° to reduce the effect of spatial autocorrelation. A total of 20 and 28 localities were used for A. bonvaloti and A. fuliginosus, respectively. Eighty per cent of the occurrence data were randomly selected to train the models and the remaining 20% were used to test the models. We set the number of maximum iterations to 2000 and the number of replicates to 10. The values of AUC and binomial probabilities were used to assess the models’ performance. To aid model interpretation, we used a conservative threshold—the lowest presence threshold[91] to distinguish ‘suitable’ from ‘unsuitable’ areas. The output maps were imported to ARC GIS 9.3 (ESRI, Redlands, CA) to produce final maps.

Contemporary gene exchange

We used a Bayesian clustering method, implemented in NEWHYBRIDS[37] to assign to each individual the posterior probabilities of being pure A. bonvaloti, pure A. fuliginosus, F1 hybrid, F2 hybrid, backcross with A. bonvaloti, or backcross with A. fuliginosus. We ran the MCMC chain for 0.4 million iterations, with the first 10% as burn-in.

References

  1. Avise JC: Molecular Markers: Natural History and Evolution. 1994, New York: Chapman and Hall

    Book  Google Scholar 

  2. Avise JC: Phylogeography: the History and Formation of Species. 2000, Cambridge: Harvard University Press

    Google Scholar 

  3. Edwards SV, Kingan SB, Calkins JD, Balakrishnan CN, Jennings WB, Swanson WJ, Sorenson MD: Speciation in birds: Genes, geography, and sexual selection. Proc Natl Acad Sci U S A. 2005, 102: 6550-6557. 10.1073/pnas.0501846102.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  4. Funk DJ, Omland KE: Species-level paraphyly and polyphyly: frequency, causes, and consequences, with insights from animal mitochondrial DNA. Annu Rev Ecol Evol Syst. 2003, 34: 397-423. 10.1146/annurev.ecolsys.34.011802.132421.

    Article  Google Scholar 

  5. Maddison WP: Gene trees in species trees. Syst Biol. 1997, 46: 523-536. 10.1093/sysbio/46.3.523.

    Article  Google Scholar 

  6. Page RD, Charleston MA: Trees within trees: phylogeny and historical associations. Trends Ecol Evol. 1998, 13: 356-359. 10.1016/S0169-5347(98)01438-4.

    Article  PubMed  CAS  Google Scholar 

  7. McKay BD, Zink RM: The causes of mitochondrial DNA gene tree paraphyly in birds. Mol Phylogenet Evol. 2010, 54: 647-650. 10.1016/j.ympev.2009.08.024.

    Article  PubMed  Google Scholar 

  8. Rheindt FE, Edwards SV: Genetic introgression: an integral but neglected component of speciation in birds. Auk. 2011, 128: 620-632. 10.1525/auk.2011.128.4.620.

    Article  Google Scholar 

  9. Arnold ML: Natural Hybridization and Evolution. 1997, Oxford: Oxford University Press

    Google Scholar 

  10. Holder MT, Anderson JA, Holloway AK: Difficulties in detecting hybridization. Syst Biol. 2001, 50: 978-982. 10.1080/106351501753462911.

    Article  PubMed  CAS  Google Scholar 

  11. Morgan K, Linton Y-M, Somboon P, Saikia P, Dev V, Socheat D, Walton C: Inter-specific gene flow dynamics during the Pleistocene-dated speciation of forest-dependent mosquitoes in Southeast Asia. Mol Ecol. 2010, 19: 2269-2285. 10.1111/j.1365-294X.2010.04635.x.

    Article  PubMed  CAS  Google Scholar 

  12. Qu Y, Zhang R, Quan Q, Song G, Li SH, Lei F: Incomplete lineage sorting or secondary admixture: disentangling historical divergence from recent gene flow in the Vinous-throated parrotbill (Paradoxornis webbianus). Mol Ecol. 2012, 21: 6117-6133. 10.1111/mec.12080.

    Article  PubMed  Google Scholar 

  13. Irwin DE, Rubtsov AS, Panov EN: Mitochondrial introgression and replacement between yellowhammers (Emberiza citrinella) and pine buntings (Emberiza leucocephalos) (Aves: Passeriformes). Biol J Linn Soc. 2009, 98: 422-438. 10.1111/j.1095-8312.2009.01282.x.

    Article  Google Scholar 

  14. Melo-Ferreira J, Boursot P, Suchentrunk F, Ferrand N, Alves PC: Invasion from the cold past: extensive introgression of mountain hare (Lepus timidus) mitochondrial DNA into three other hare species in northern Iberia. Mol Ecol. 2005, 14: 2459-2464. 10.1111/j.1365-294X.2005.02599.x.

    Article  PubMed  CAS  Google Scholar 

  15. Peters JL, Zhuravlev Y, Fefelov I, Logie A, Omland KE: Nuclear loci and coalescent methods support ancient hybridization as cause of mitochondrial paraphyly between gadwall and falcated duck (Anas spp.). Evolution. 2007, 61: 1992-2006. 10.1111/j.1558-5646.2007.00149.x.

    Article  PubMed  CAS  Google Scholar 

  16. Weckstein JD, Zink RM, Blackwell-Rago RC, Nelson DA: Anomalous variation in mitochondrial genomes of White-crowned (Zonotrichia leucophrys) and Golden-crowned (Z. atricapilla) sparrows: pseudogenes, hybridization, or incomplete lineage sorting?. Auk. 2001, 118: 231-236. 10.1642/0004-8038(2001)118[0231:AVIMGO]2.0.CO;2.

    Article  Google Scholar 

  17. Alström P, Olsson U, Lei F, Wang H-t, Gao W, Sundberg P: Phylogeny and classification of the Old World Emberizini (Aves, Passeriformes). Mol Phylogenet Evol. 2008, 47: 960-973. 10.1016/j.ympev.2007.12.007.

    Article  PubMed  Google Scholar 

  18. Maddison WP, Knowles LL: Inferring phylogeny despite incomplete lineage sorting. Syst Biol. 2006, 55: 21-30. 10.1080/10635150500354928.

    Article  PubMed  Google Scholar 

  19. Ballard JWO, Whitlock MC: The incomplete natural history of mitochondria. Mol Ecol. 2004, 13: 729-744. 10.1046/j.1365-294X.2003.02063.x.

    Article  PubMed  Google Scholar 

  20. Carling MD, Brumfield RT: Haldane’s rule in an avian system: using cline theory and divergence population genetics to test for differential introgression of mitochondrial, autosomal, and sex-linked loci across the Passerina bunting hybrid zone. Evolution. 2008, 62: 2600-2615. 10.1111/j.1558-5646.2008.00477.x.

    Article  PubMed  CAS  Google Scholar 

  21. Hey J: On the number of New World founders: a population genetic portrait of the peopling of the Americas. PLoS Biol. 2005, 3: e193-10.1371/journal.pbio.0030193.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Hey J, Nielsen R: Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis. Genetics. 2004, 167: 747-760. 10.1534/genetics.103.024182.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  23. Nielsen R, Wakeley J: Distinguishing migration from isolation: a Markov chain Monte Carlo approach. Genetics. 2001, 158: 885-896.

    PubMed  CAS  PubMed Central  Google Scholar 

  24. Backström N, Väli Ü: Sex- and species-biased gene flow in a spotted eagle hybrid zone. BMC Evol Biol. 2011, 11: 100-10.1186/1471-2148-11-100.

    Article  PubMed  PubMed Central  Google Scholar 

  25. Melo-Ferreira J, Boursot P, Carneiro M, Esteves PJ, Farelo L, Alves PC: Recurrent introgression of mitochondrial DNA among hares (Lepus spp.) revealed by species-tree inference and coalescent simulations. Syst Biol. 2011, 61: 367-381.

    Article  PubMed  Google Scholar 

  26. Warren BH, Bermingham E, Bourgeois Y, Estep LK, Prys-Jones RP, Strasberg D, Thébaud C: Hybridization and barriers to gene flow in an island bird radiation. Evolution. 2012, 66: 1490-1505. 10.1111/j.1558-5646.2011.01550.x.

    Article  PubMed  Google Scholar 

  27. Carling MD, Brumfield RT: Speciation in Passerina buntings: introgression patterns of sex-linked loci identify a candidate gene region for reproductive isolation. Mol Ecol. 2009, 18: 834-847. 10.1111/j.1365-294X.2008.04038.x.

    Article  PubMed  CAS  Google Scholar 

  28. Storchová R, Reif J, Nachman MW: Female heterogamety and speciation: reduced introgression of the Z chromosome between two species of nightingales. Evolution. 2009, 64: 456-471.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Dai C, Chen K, Zhang R, Yang X, Yin Z, Tian H, Zhang Z, Hu Y, Lei F: Molecular phylogenetic analysis among species of Paridae, Remizidae and Aegithalos based on mtDNA sequences of COI and cyt b. Chinese Birds. 2010, 1: 112-123. 10.5122/cbirds.2010.0003.

    Article  Google Scholar 

  30. Päckert M, Martens J, Sun Y-H: Phylogeny of long-tailed tits and allies inferred from mitochondrial and nuclear markers (Aves: Passeriformes, Aegithalidae). Mol Phylogenet Evol. 2010, 55: 952-967. 10.1016/j.ympev.2010.01.024.

    Article  PubMed  Google Scholar 

  31. del Hoyo J, Elliott A, Christie D: Handbook of the Birds of the World. Penduline Tits to Shrikes. 2008, Barcelona: Lynx Edicions

    Google Scholar 

  32. Harrap S, Quinn D: Chickadees, tits, nuthatches & treecreepers. 1996, Princeton: Princeton University Press

    Google Scholar 

  33. Li G, Zheng B, Liu G: Passeriformes: Paridae - Zosteropidae. Fauna Sinica. 1982, Beijing: Science Press

    Google Scholar 

  34. Sorenso MD, Quinn TW: Numts: a challenge for avian systematics and population biology. Auk. 1998, 115: 214-221. 10.2307/4089130.

    Article  Google Scholar 

  35. Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123: 585-595.

    PubMed  CAS  PubMed Central  Google Scholar 

  36. Fu YX, Li WH: Statistical tests of neutrality of mutations. Genetics. 1993, 133: 693-709.

    PubMed  CAS  PubMed Central  Google Scholar 

  37. Anderson EC, Thompson EA: A model-based method for identifying species hybrids using multilocus genetic data. Genetics. 2002, 160: 1217-1229.

    PubMed  CAS  PubMed Central  Google Scholar 

  38. Axelsson E, Smith NGC, Sundström H, Berlin S, Ellegren H: Male-biased mutation rate and divergence in autosomal, Z-linked and W-linked introns of chicken and turkey. Mol Biol Evol. 2004, 21: 1538-1547. 10.1093/molbev/msh157.

    Article  PubMed  CAS  Google Scholar 

  39. Weir JT, Schluter D: Calibrating the avian molecular clock. Mol Ecol. 2008, 17: 2321-2328. 10.1111/j.1365-294X.2008.03742.x.

    Article  PubMed  CAS  Google Scholar 

  40. Zink RM, Barrowclough GF: Mitochondrial DNA under siege in avian phylogeography. Mol Ecol. 2008, 17: 2107-2121. 10.1111/j.1365-294X.2008.03737.x.

    Article  PubMed  CAS  Google Scholar 

  41. Morando M, Avila LJ, Baker J, Sites JW: Phylogeny and phylogeography of the Liolaemus darwinii complex (Squamata: Liolaemidae): evidence for introgression and incomplete lineage sorting. Evolution. 2004, 58: 842-861. 10.1111/j.0014-3820.2004.tb00416.x.

    Article  PubMed  Google Scholar 

  42. Ng J, Glor RE: Genetic differentiation among populations of a Hispaniolan trunk anole that exhibit geographical variation in dewlap colour. Mol Ecol. 2011, 20: 4302-4317. 10.1111/j.1365-294X.2011.05267.x.

    Article  PubMed  Google Scholar 

  43. Wirtz P: Mother species–father species: unidirectional hybridization in animals with female choice. Anim Behav. 1999, 58: 1-12. 10.1006/anbe.1999.1144.

    Article  PubMed  Google Scholar 

  44. Stein AC, Uy JAC: Unidirectional introgression of a sexually selected trait across an avian hybrid zone: a role for female choice?. Evolution. 2006, 60: 1476-1485.

    PubMed  Google Scholar 

  45. Ryan MJ, Wagner WE: Asymmetries in mating preferences between species: female swordtails prefer heterospecific males. Science. 1987, 236: 595-597. 10.1126/science.236.4801.595.

    Article  PubMed  CAS  Google Scholar 

  46. Krosby M, Rohwer S: A 2000 km genetic wake yields evidence for northern glacial refugia and hybrid zone movement in a pair of songbirds. Proc Biol Sci. 2009, 276: 615-621. 10.1098/rspb.2008.1310.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  47. McDonald DB, Clay RP, Brumfield RT, Braun MJ: Sexual selection on plumage and behavior in an avian hybrid zone: experimental tests of male‒male interactions. Evolution. 2001, 55: 1443-1451. 10.1111/j.0014-3820.2001.tb00664.x.

    Article  PubMed  CAS  Google Scholar 

  48. Clarke AL, Sæther BE, Røskaft E: Sex biases in avian dispersal: a reappraisal. Oikos. 1997, 79: 429-438. 10.2307/3546885.

    Article  Google Scholar 

  49. Greenwood PJ: Mating systems, philopatry and dispersal in birds and mammals. Anim Behav. 1980, 28: 1140-1162. 10.1016/S0003-3472(80)80103-5.

    Article  Google Scholar 

  50. Good JM, Hird S, Reid N, Demboski JR, Steppan SJ, Martin-Nims TR, Sullivan J: Ancient hybridization and mitochondrial capture between two species of chipmunks. Mol Ecol. 2008, 17: 1313-1327. 10.1111/j.1365-294X.2007.03640.x.

    Article  PubMed  CAS  Google Scholar 

  51. Barbanera F, Zuffi MAL, Guerrini M, Gentilli A, Tofanelli S, Fasola M, Dini F: Molecular phylogeography of the asp viper Vipera aspis (Linnaeus, 1758) in Italy: Evidence for introgressive hybridization and mitochondrial DNA capture. Mol Phylogenet Evol. 2009, 52: 103-114. 10.1016/j.ympev.2009.02.006.

    Article  PubMed  CAS  Google Scholar 

  52. Willis SC, Farias IP, Ortí G: Testing mitochondrial capture and deep coalescence in Amazonian cichlid fishes (Cichlidae: Cichla). Evolution. 2013, 68: 256-268.

    Article  PubMed  Google Scholar 

  53. Toews DPL, Brelsford A: The biogeography of mitochondrial and nuclear discordance in animals. Mol Ecol. 2012, 21: 3907-3930. 10.1111/j.1365-294X.2012.05664.x.

    Article  PubMed  CAS  Google Scholar 

  54. Chan KMA, Levin SA: Leaky prezygotic isolation and porous genomes: rapid introgression of maternally inherited DNA. Evolution. 2005, 59: 720-729. 10.1111/j.0014-3820.2005.tb01748.x.

    Article  PubMed  CAS  Google Scholar 

  55. Harrison RG: Hybrid zones: windows on evolutionary process. Oxf Surv Evol Biol. 1990, 7: 69-128.

    Google Scholar 

  56. Arnold J: Cytonuclear disequilibria in hybrid zones. Annu Rev Ecol Evol Syst. 1993, 24: 521-554. 10.1146/annurev.es.24.110193.002513.

    Article  Google Scholar 

  57. Doiron S, Bernatchez L, Blier PU: A comparative mitogenomic analysis of the potential adaptive value of arctic charr mtDNA introgression in brook charr populations (Salvelinus fontinalis Mitchill). Mol Phylogenet Evol. 2002, 19: 1902-1909.

    Article  CAS  Google Scholar 

  58. Cheviron ZA, Brumfield RT: Migration-selection balance and local adaptation of mitochondrial haplotypes in Rufous-collared Sparrows (Zonotrichia Capensis) along an elevational gradient. Evolution. 2009, 63: 1593-1605. 10.1111/j.1558-5646.2009.00644.x.

    Article  PubMed  Google Scholar 

  59. Ehinger M, Fontanillas P, Petit E, Perrin N: Mitochondrial DNA variation along an altitudinal gradient in the greater white‒toothed shrew, Crocidura russula. Mol Ecol. 2002, 11: 939-945. 10.1046/j.1365-294X.2002.01487.x.

    Article  PubMed  CAS  Google Scholar 

  60. Wise CA, Sraml M, Easteal S: Departure from neutrality at the mitochondrial NADH dehydrogenase subunit 2 gene in humans, but not in chimpanzees. Genetics. 1998, 148: 409-421.

    PubMed  CAS  PubMed Central  Google Scholar 

  61. Zink RM: Natural selection on mitochondrial DNA in Parus and its relevance for phylogeographic studies. Proc Biol Sci. 2005, 272: 71-78. 10.1098/rspb.2004.2908.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  62. Ferris SD, Sage RD, Huang C-M, Nielsen JT, Ritte U, Wilson AC: Flow of mitochondrial DNA across a species boundary. Proc Natl Acad Sci U S A. 1983, 80: 2290-2294. 10.1073/pnas.80.8.2290.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  63. Tang Q-Y, Liu S-Q, Yu D, Liu H-Z, Danley PD: Mitochondrial capture and incomplete lineage sorting in the diversification of balitorine loaches (Cypriniformes, Balitoridae) revealed by mitochondrial and nuclear genes. Zool Scr. 2012, 41: 233-247. 10.1111/j.1463-6409.2011.00530.x.

    Article  Google Scholar 

  64. Waltari E, Guralnick RP: Ecological niche modelling of montane mammals in the Great Basin, North America: examining past and present connectivity of species across basins and ranges. J Biogeography. 2009, 36: 148-161. 10.1111/j.1365-2699.2008.01959.x.

    Article  Google Scholar 

  65. Wang W, McKay BD, Dai C, Zhao N, Zhang R, Qu Y, Song G, Li S-H, Liang W, Yang X, Pasquet E, Lei F: Glacial expansion and diversification of an East Asian montane bird, the green-backed tit (Parus monticolus). J Biogeography. 2013, 40: 1156-1169. 10.1111/jbi.12055.

    Article  Google Scholar 

  66. Galbreath KE, Hafner DJ, Zamudio KR, Agnew K: Isolation and introgression in the Intermountain West: contrasting gene genealogies reveal the complex biogeographic history of the American pika (Ochotona princeps). J Biogeography. 2010, 37: 344-362. 10.1111/j.1365-2699.2009.02201.x.

    Article  Google Scholar 

  67. Zieliński P, Nadachowska-Brzyska K, Wielstra B, Szkotak R, Covaciu-Marcov SD, Cogălniceanu D, Babik W: No evidence for nuclear introgression despite complete mtDNA replacement in the Carpathian newt (Lissotriton montandoni). Mol Ecol. 2013, 22: 1884-1903. 10.1111/mec.12225.

    Article  PubMed  Google Scholar 

  68. Birckhead H, Sage D, Sheldon WG, Carter TD: The Birds of the Sage West China Expedition. 1937, New York: The American Museum of Natural History

    Google Scholar 

  69. Sorenson MD, Ast JC, Dimcheff DE, Yuri T, Mindell DP: Primers for a PCR-based approach to mitochondrial genome sequencing in birds and other vertebrates. Mol Phylogenet Evol. 1999, 12: 105-114. 10.1006/mpev.1998.0602.

    Article  PubMed  CAS  Google Scholar 

  70. Kimball RT, Braun EL, Barker FK, Bowie RCK, Braun MJ, Chojnowski JL, Hackett SJ, Han K-L, Harshman J, Heimer-Torres V, Holznagel W, Huddleston CJ, Marks BD, Miglia KJ, Moore WS, Reddy S, Sheldon FH, Smith JV, Witt CC, Yuri T: A well-tested set of primers to amplify regions spread across the avian genome. Mol Phylogenet Evol. 2009, 50: 654-660. 10.1016/j.ympev.2008.11.018.

    Article  PubMed  CAS  Google Scholar 

  71. Irestedt M, Johansson US, Parsons TJ, Ericson PGP: Phylogeny of major lineages of suboscines (Passeriformes) analysed by nuclear DNA sequence data. J Avian Biol. 2001, 32: 15-25. 10.1034/j.1600-048X.2001.320103.x.

    Article  Google Scholar 

  72. Borge T, Webster MT, Andersson G, Saetre GP: Contrasting patterns of polymorphism and divergence on the Z chromosome and autosomes in two Ficedula flycatcher species. Genetics. 2005, 171: 1861-1873. 10.1534/genetics.105.045120.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  73. Fridolfsson AK, Ellegren H: A simple and universal method for molecular sexing of non-ratite birds. J Avian Biol. 1999, 30: 116-121. 10.2307/3677252.

    Article  Google Scholar 

  74. Tamura K, Dudley J, Nei M, Kumar S: MEGA4: molecular evolutionary genetics analysis (MEGA) software version 4.0. Mol Biol Evol. 2007, 24: 1596-1599. 10.1093/molbev/msm092.

    Article  PubMed  CAS  Google Scholar 

  75. Stephens M, Smith NJ, Donnelly P: A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001, 68: 978-989. 10.1086/319501.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  76. Carling MD, Lovette IJ, Brumfield RT: Historical divergence and gene flow: coalescent analyses of mitochondrial, autosomal and sex-linked loci in Passerina buntings. Evolution. 2010, 64: 1762-1772. 10.1111/j.1558-5646.2010.00954.x.

    Article  PubMed  Google Scholar 

  77. Harrigan RJ, Mazza ME, Sorenson MD: Computation vs. cloning: evaluation of two methods for haplotype determination. Mol Ecol Resour. 2008, 8: 1239-1248. 10.1111/j.1755-0998.2008.02241.x.

    Article  PubMed  CAS  Google Scholar 

  78. Librado P, Rozas J: DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009, 25: 1451-1452. 10.1093/bioinformatics/btp187.

    Article  PubMed  CAS  Google Scholar 

  79. Hey J: HKA. [https://bio.cst.temple.edu/~hey/software/software.htm]

  80. Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000, 9: 1657-1659. 10.1046/j.1365-294x.2000.01020.x.

    Article  PubMed  CAS  Google Scholar 

  81. Drummond AJ, Rambaut A: BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007, 7: 214-10.1186/1471-2148-7-214.

    Article  PubMed  PubMed Central  Google Scholar 

  82. Nylander J: MrModeltest v2. Program Distributed by the Author. 2004, Uppsala University: Evolutionary Biology Centre

    Google Scholar 

  83. Hasegawa M, Kishino H, Yano T-a: Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 1985, 22: 160-174. 10.1007/BF02101694.

    Article  PubMed  CAS  Google Scholar 

  84. Tracer V1.5. [http://tree.bio.ed.ac.uk/software/tracer/]

  85. FigTree version 1.3. 1. [http://tree.bio.ed.ac.uk/software/figtree/]

  86. Hey J, Nielsen R: Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proc Natl Acad Sci U S A. 2007, 104: 2785-2790. 10.1073/pnas.0611164104.

    Article  PubMed  CAS  PubMed Central  Google Scholar 

  87. Strasburg JL, Rieseberg LH: How robust are “Isolation with Migration” analyses to violations of the IM model? A simulation study. Mol Biol Evol. 2009, 27: 297-310.

    Article  PubMed  PubMed Central  Google Scholar 

  88. Excoffier L, Laval G, Schneider S: Arlequin (version 3.0): an integrated software package for population genetics data analysis. Evol Bioinformatics Online. 2005, 1: 47-50.

    CAS  Google Scholar 

  89. Phillips SJ, Anderson RP, Schapire RE: Maximum entropy modelling of species geographic distributions. Ecol Model. 2006, 190: 231-259. 10.1016/j.ecolmodel.2005.03.026.

    Article  Google Scholar 

  90. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A: Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005, 25: 1965-1978. 10.1002/joc.1276.

    Article  Google Scholar 

  91. Pearson RG, Raxworthy CJ, Nakamura M, Peterson TA: Predicting species distributions from small numbers of occurrence records: a test case using cryptic geckos in Madagascar. J Biogeography. 2007, 34: 102-117.

    Article  Google Scholar 

Download references

Acknowledgements

We thank Zuohua Yin, Huatao Liu and Bin Gao for field assistance and sample collection, Yongjie Wu, Baoyan Liu and Ruiying Zhang for data analyses assistance, and Jiakuan Chen, Jun Rong, Shane DuBay, Bailey D. McKay, Janette Norman and two anonymous reviewers for providing valuable comments. We also thank the National Zoological Museum of China for providing sampling site information of the two species. This work was supported by the State Key Program of National Natural Science of China (No 31330073), Major International (Regional) Joint Research Project (No: 31010103901) to FL and PA, Chinese Academy of Sciences Visiting Professorship for Senior International Scientists (No. 2011T2S04) to PA, CAS-IOZ Innovation Program (KSCX2-EW-J-2) to FL and by NSFC program J1210002 to WW.

Author information

Authors and Affiliations

  1. Key Laboratory of Zoological Systematics and Evolution, Institute of Zoology, Chinese Academy of Sciences, No.1 Beichen West Road, Chaoyang District, 100101, Beijing, China

    Wenjuan Wang, Per Alström, Chunlan Zhang, Yanhua Qu, Gang Song & Fumin Lei

  2. Center for Watershed Ecology, Institute of Life Science, Nanchang University, 330031, Nanchang, China

    Wenjuan Wang

  3. School of Chemistry and Life Sciences, Guizhou Normal College, 550001, Guiyang, China

    Chuanyin Dai

  4. Swedish Species Information Centre, Swedish University of Agricultural Sciences, SE-75007, Uppsala, Sweden

    Per Alström

  5. Department of Life Science, National Taiwan Normal University, Taipei, Taiwan

    Shou-Hsien Li

  6. State Key Laboratory of Genetic Resources and Evolution, Kunming Institute of Zoology, Chinese Academy of Sciences, 650223, Kunming, China

    Xiaojun Yang

  7. Key Laboratory of Ecological Impacts of Hydraulic-Projects and Restoration of Aquatic Ecosystem of Ministry of Water Resources, Institute of Hydroecology, Ministry of Water Resources and Chinese Academy of Sciences, 430079, Wuhan, China

    Na Zhao

Authors
  1. Wenjuan Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Chuanyin Dai
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Per Alström
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Chunlan Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Yanhua Qu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Shou-Hsien Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Xiaojun Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Na Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Gang Song
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Fumin Lei
    View author publications

    You can also search for this author in PubMed Google Scholar

Corresponding author

Correspondence to Fumin Lei.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

FL, CD and WW conceived and designed the experiments; CD, XY and GS collected samples; WW, CD and NZ performed the experiments; WW, CD and CZ analyzed the data; WW, PA and FL wrote the paper; YQ and SL provided valuable suggestions. All authors read and approved the final manuscript.

Electronic supplementary material

Authors’ original submitted files for images

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver ( https://creativecommons.org/publicdomain/zero/1.0/ ) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, W., Dai, C., Alström, P. et al. Past hybridization between two East Asian long-tailed tits (Aegithalos bonvaloti and A. fuliginosus). Front Zool 11, 40 (2014). https://0-doi-org.brum.beds.ac.uk/10.1186/1742-9994-11-40

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://0-doi-org.brum.beds.ac.uk/10.1186/1742-9994-11-40

Keywords

Frontiers in Zoology

ISSN: 1742-9994