- Open Access
Past hybridization between two East Asian long-tailed tits (Aegithalos bonvaloti and A. fuliginosus)
Frontiers in Zoology volume 11, Article number: 40 (2014)
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.
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.
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.
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 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. 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. Massive hybridization over a long period could also homogenize nuDNA, making conclusions more difficult to draw. 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. 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. 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. suggested that A. bonvaloti and A. fuliginosus may “represent just distinct morphotypes at an intraspecific differentiation level”. However, both Dai et al. and Päckert et al. 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.
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.
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).
The absence of ambiguous sites and stop codons suggested that these mtDNA sequences were not “numts”. 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 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.
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.
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 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.
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.
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).
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.
Contemporary gene exchange
We used NEWHYBRIDS 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).
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, display more attractive courtship, 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. 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, and has been suggested to be the source of shallow genetic structure in Willow Tit (Parus montanus;). 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.
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.
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. 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, R50/R51, BRM15F/BRM15R, F2550/R2718, and VLDLR9F/VLDLR9R. 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. 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. 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 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. DnaSP was used to calculate polymorphic sites, number of haplotypes, haplotype diversities and nucleotide diversities.
We used TCS 1.21 to produce a maximum-parsimony network for each locus. The Bayesian algorithm implemented in BEAST 1.7.0 was used to reconstruct phylogenetic relationships among mtDNA haplotypes. We selected the best-fit model of nucleotide substitution using the program MrModeltest 2.3. According to the Akaike information criterion (AIC), the best fit model for mtDNA was the HKY model. 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. The output trees file was summarized into a maximum clade credibility tree using TreeAnnotator of the BEAST package. FigTree 1.3.1 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.
In order to distinguish between ILS and hybridization, we applied the coalescent-based IM model implemented in the program IMa. 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. 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.
We conducted mismatch distribution analysis using Arlequin 3.5. 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 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 (; 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 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 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.
Avise JC: Molecular Markers: Natural History and Evolution. 1994, New York: Chapman and Hall
Avise JC: Phylogeography: the History and Formation of Species. 2000, Cambridge: Harvard University Press
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.
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.
Maddison WP: Gene trees in species trees. Syst Biol. 1997, 46: 523-536. 10.1093/sysbio/46.3.523.
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.
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.
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.
Arnold ML: Natural Hybridization and Evolution. 1997, Oxford: Oxford University Press
Holder MT, Anderson JA, Holloway AK: Difficulties in detecting hybridization. Syst Biol. 2001, 50: 978-982. 10.1080/106351501753462911.
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.
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.
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.
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.
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.
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.
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.
Maddison WP, Knowles LL: Inferring phylogeny despite incomplete lineage sorting. Syst Biol. 2006, 55: 21-30. 10.1080/10635150500354928.
Ballard JWO, Whitlock MC: The incomplete natural history of mitochondria. Mol Ecol. 2004, 13: 729-744. 10.1046/j.1365-294X.2003.02063.x.
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.
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.
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.
Nielsen R, Wakeley J: Distinguishing migration from isolation: a Markov chain Monte Carlo approach. Genetics. 2001, 158: 885-896.
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.
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.
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.
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.
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.
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.
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.
del Hoyo J, Elliott A, Christie D: Handbook of the Birds of the World. Penduline Tits to Shrikes. 2008, Barcelona: Lynx Edicions
Harrap S, Quinn D: Chickadees, tits, nuthatches & treecreepers. 1996, Princeton: Princeton University Press
Li G, Zheng B, Liu G: Passeriformes: Paridae - Zosteropidae. Fauna Sinica. 1982, Beijing: Science Press
Sorenso MD, Quinn TW: Numts: a challenge for avian systematics and population biology. Auk. 1998, 115: 214-221. 10.2307/4089130.
Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123: 585-595.
Fu YX, Li WH: Statistical tests of neutrality of mutations. Genetics. 1993, 133: 693-709.
Anderson EC, Thompson EA: A model-based method for identifying species hybrids using multilocus genetic data. Genetics. 2002, 160: 1217-1229.
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.
Weir JT, Schluter D: Calibrating the avian molecular clock. Mol Ecol. 2008, 17: 2321-2328. 10.1111/j.1365-294X.2008.03742.x.
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.
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.
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.
Wirtz P: Mother species–father species: unidirectional hybridization in animals with female choice. Anim Behav. 1999, 58: 1-12. 10.1006/anbe.1999.1144.
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.
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.
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.
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.
Clarke AL, Sæther BE, Røskaft E: Sex biases in avian dispersal: a reappraisal. Oikos. 1997, 79: 429-438. 10.2307/3546885.
Greenwood PJ: Mating systems, philopatry and dispersal in birds and mammals. Anim Behav. 1980, 28: 1140-1162. 10.1016/S0003-3472(80)80103-5.
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.
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.
Willis SC, Farias IP, Ortí G: Testing mitochondrial capture and deep coalescence in Amazonian cichlid fishes (Cichlidae: Cichla). Evolution. 2013, 68: 256-268.
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.
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.
Harrison RG: Hybrid zones: windows on evolutionary process. Oxf Surv Evol Biol. 1990, 7: 69-128.
Arnold J: Cytonuclear disequilibria in hybrid zones. Annu Rev Ecol Evol Syst. 1993, 24: 521-554. 10.1146/annurev.es.24.110193.002513.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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
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.
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.
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.
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.
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.
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.
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.
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.
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.
Librado P, Rozas J: DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009, 25: 1451-1452. 10.1093/bioinformatics/btp187.
Hey J: HKA. [https://bio.cst.temple.edu/~hey/software/software.htm]
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.
Drummond AJ, Rambaut A: BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol. 2007, 7: 214-10.1186/1471-2148-7-214.
Nylander J: MrModeltest v2. Program Distributed by the Author. 2004, Uppsala University: Evolutionary Biology Centre
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.
Tracer V1.5. [http://tree.bio.ed.ac.uk/software/tracer/]
FigTree version 1.3. 1. [http://tree.bio.ed.ac.uk/software/figtree/]
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.
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.
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.
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.
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.
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.
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.
The authors declare that they have no competing interests.
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
Below are the links to the authors’ original submitted files for images.
About this article
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://doi.org/10.1186/1742-9994-11-40
- Incomplete lineage sorting
- Multilocus analyses
- Ecological niche modelling