Introgressive hybridization among Chinese Lepus species
As seen from the analyses of mtDNA and nuclear DNA sequences of Chinese hares, most of the phylogenetic conflicts between the mtDNA and nuclear DNA gene trees suggest multiple episodes of introgression. Alternatively, incomplete lineage sorting of the nuclear gene may also result in the observed "mixing" of the nuclear gene haplotypes; however, most of the samples with the discrepancy between the mitochondrial and nuclear genotypes have congruence between their morphology and nuclear genotype except for some L. capensis and L. sinensis individuals (Additional file 2). Although some of the L. capensis and L. sinensis individuals have incongruence between their morphology and nuclear genotype, the "L. capensis group" and the "L. sinensis group" is divergent enough to distinguish each other (Figure 2), therefore, the nuclear MGF gene is informative enough to distinguish Chinese hare species and the observed "mixing" of nuclear gene haplotypes likely is due to introgressive hybridization rather than incomplete lineage sorting. In addition, the observation that a large number of the heterozygous genotypes have alleles from divergent clades (Additional file 2) is also inconsistent with the alternative hypothesis of incomplete lineage sorting. The evidence suggests that introgression has occurred between seven pairs of the hare species examined here.
L. sinensis × L. capensis, L. sinensis (♀) × L. mandshuricus (♂), L. sinensis × L. timidus
Five L. capensis individuals (CA9, CA11- CA13 and CA36) were placed within the L. sinensis lineage on the mtDNA tree (Figure 1). These results likely reflect recent introgression of L. sinensis mtDNA into L. capensis; however, results from the nuclear gene analysis suggested bidirectional and recent introgression between these species, as evidenced by the observation that the MGF haplotypes from L. capensis and L. sinensis were admixed and fell within two distantly related groups. Furthermore, one haplotype (Hap3) that nested in the L. capensis group and two haplotypes (Hap49, 53) that nested in the L. sinensis group were shared between L. capensis and L. sinensis (Figure 2). In addition, ten L. capensis individuals (CA1, CA2, CA10, CA17-CA19, CA21, CA22, CA29 and CA31) and two L. sinensis individuals (S2 and S3) were identified as heterozygotes at the MGF locus for two alleles from L. capensis and L. sinensis (Additional file 2). The available data, however, cannot eliminate the possibility that some of the heterozygotes were produced by ancient hybridization events.
Five individuals (M8, M10, M12, M13 and M18) of L. mandshuricus were grouped into the L. sinensis lineage on the mtDNA tree (Figure 1), while their MGF haplotypes fell within the L. mandshuricus clade (Figure 2). This can be explained by a unidirectional mtDNA introgression from L. sinensis into L. mandshuricus. In addition, two individuals of L. mandshuricus (M10 and M12) have low sequence divergence relative to the L. sinensis individuals, thus supporting a recent age for the hybridization events. The divergence time estimation indicated that the remaining three L. mandshuricus individuals (M8, M13 and MA18) shared a most recent ancestor with the L. sinensis samples at ca. 0.414 MYA (± 0.133 MYA) (Figure 1). This would suggest that the mitochondrial genome of L. sinensis introgressed into L. mandshuricus no later than 0.414MYA. Our estimate of the time of introgression may be biased somewhat if additional extant lineages were not sampled, or if there have been extinction events. Indeed, such a bias has been inferred in the case of Lepidiolamprologus . However, in the present study three ancient hybridization events, involving five species (Figure 1; also see below), were detected. This strongly supports the hypothesis of ancient hybridization and not solely recent introgression events.
Two individuals of L. timidus, T10 and T11, nested within the L. sinensis clade on the mtDNA tree (Figure 1), while their MGF haplotypes grouped within the L. timidus clade (Figure 2), once again suggesting mtDNA introgression (in this case, from L. sinensis into L. timidus). In addition, the detection of one L. timidus haplotype (Hap65) within the L. sinensis group on the MGF gene tree (Figure 2) was consistent with the occurrence of nuclear DNA introgression from L. timidus into L. sinensis. Taken together, the mtDNA and nuclear DNA analyses supported reciprocal introgression between L. timidus and L. sinensis. An alternative explanation would be that the introgressive hybridization of L. sinensis × L. capensis, L. sinensis (♀) × L. mandshuricus (♂), and L. sinensis (♀) ×L. timidus (♂) resulted from an initial hybridization between L. sinensis and L. capensis, L. mandshuricus or L. timidus, followed by a second hybridization among L. capensis, L. mandshuricus and L. timidus. Supporting this latter hypothesis was the inference of hybridization between L. timidus and L. mandshuricus as well as L. timidus and L. capensis (see below).
L. timidus (♀) × L. mandshuricus (♂), L. timidus (♀) × L. capensis (♂)
Alves et al.  suggested that mtDNA introgression from L. timidus might occur in several other regions in addition to Iberia and our results confirm this hypothesis. In the L. timidus mtDNA lineage, nine L. mandshuricus individuals (M1-M4, M7, M9, M11, M14 and M19) and three L. capensis individuals (CA14, CA15 and CA35) were included (Figure 1). The MGF haplotypes of the nine L. mandshuricus individuals were grouped with the other L. mandshuricus haplotypes except for Hap38, which is nested within the L. timidus clade. The nuclear sequences of the three L. capensis individuals were placed in either the L. capensis group (CA15) or the L. sinensis group (CA35) (Additional file 2). In addition, one L. timidus haplotype (Hap7) was placed into the L. capensis group and two haplotypes (Hap15 and Hap21) were shared with L. capensis on the MGF gene tree (Figure 2). Our results thus suggest that the L. mandshuricus and L. capensis individuals likely reflect the capture of L. timidus mitochondrial genomes via several hybridization events between females of L. timidus and males of L. mandshuricus or L. capensis. Furthermore, the nuclear DNA analysis indicated that the hybridization events between L. timidus and L. mandshuricus or L. capensis were recent and a second hybridization event likely occurred between L. capensis and L. sinensis as well. An additional nine Cyt b sequences of L. mandshuricus accessed from Genbank (accession numbers: AY650894-AY650899, AJ279423, DQ793162, DQ793163) were also grouped into the L. timidus lineage on our Cyt b tree (data not show), further supporting the inference of mtDNA introgression from L. timidus into L. mandshuricus. Finally, one L. timidus Cyt b sequence (accession number: AY745108) collected from GenBank was grouped within the L. capensis lineage in the Cyt b phylogeny (data not show), suggesting the introgression of mtDNA from L. capensis into L. timidus. The L. capensis individual, CA35, shared a most recent common ancestor with L. timidus samples at 0.976 MYA (± 0.284 MYA) (Figure 1), indicating that the hybridization between L. timidus and L. capensis took place at about this point in time.
L. yarkandensis (♀) × L. capensis(♂), L. capensis × L. oiostolus
Several observations suggest the role of recent hybridization between L. capensis and L. yarkandensis. These observations include: one MGF haplotype (Hap4) of L. yarkandensis nested within the L. capensis group; one MGF haplotype (Hap13) of L. yarkandensis shared with L. capensis-2 (Figure 2); and two MGF alleles of four L. yarkandensis individuals (Y7, Y10, Y14, Y15) were grouped into either the L. yarkandensis and L. capensis groups (Additional file 2).
The newly identified L. capensis-2 lineage might have originated through ancient hybridization between female L. yarkandensis and male L. capensis at ca 0.988 MYA (± 0.285 MYA) (Figure 1; also see below). Interestingly, one L. capensis haplotype (Hap31) was grouped into the L. oiostolus clade (Figure 2). This may reflect a recent hybridization event between L. capensis and L. oiostolus. In addition, hybridization between L. oiostolus and L. capensis was also evidenced by the fact that one L. oiostolus Cyt b sequence (accession number: AJ279427) from GenBank was nested within the L. capensis lineage (data not shown).
The possible forces driving the introgressive hybridization among Chinese Lepus species
Compared with introgression reported previously for other species complexes [[18, 21–23, 54, 55]], especially within other Lepus species [[17, 25–27, 31, 32, 50]], Chinese hares demonstrated a previously unreported level of complexity. First, most of the previous studies within the genus Lepus reported only the presence of mtDNA introgression of L. timidus, in which unidirectional mtDNA introgression events from L. timidus into populations of L. europaeus in Sweden [[25–27]] and into L. granatensis and L. europaeus in Iberian Peninsula [[17, 31, 32]] were inferred, with the identification of nuclear introgression lacking or infrequent [[48, 50, 53]]. Similar scenarios have also been posited for several other species clades [[18, 22, 23, 54–56]]. Second, although contemporaneous introgression has been reported for several non-hare species [[21, 57, 58]], most of these cases of introgression involve at the most two or three species. In comparison, the present study demonstrated that introgression among Chinese Lepus species has most likely been a continuous and recent process involving multiple waves of hybridization among six species.
Our results also suggest that introgression has occurred in multiple directions, involving both mitochondrial and nuclear DNA. The mtDNA of L. mandshuricus was replaced entirely by that of L. timidus and L. sinensis. Although several previous studies have provided convincing evidence for complete local or even range-wide replacement of mtDNA following hybridization in insects, fishes and reptilia [[18, 23, 55, 59]], such observations are infrequent for mammals. The hybridization events among Chinese hares are possibly driven by species abundance asymmetries, mating preferences, range expansion and low divergence between Lepus species.
In the studies of Thulin & Tegelström , Hubbs  and Wirtz , they suggested that the females of a rare species (e.g. L. timidus) would generally hybridize with the males of a common species (e.g. L. europaeus). In the present study, the hybridization between females of L. timidus (rare species) and males of L. capensis (common species), and between females of L. timidus (rare species) and males of L. mandshuricus (common species) are congruent with their conclusion. However, the hybridization between females of L. sinensis (common species) and males of L. capensis (common species) is inconsistent with this hypothesis. In addition, Thulin & Tegelström  and Grant & Grant  indicated that unidirectional hybridization would usually occur between the females of smaller species (e.g. L. timidus) and the males of larger-bodied species (e.g. L. europaeus). Among Chinese hares, hybridization between females of L. sinensis (smaller species) and males of L. capensis (larger species), L. timidus (larger species) or L. mandshuricus (larger species) support this hypothesis. In contrast, the hybridizations between female L. timidus (larger species) and male L. capensis (smaller species) or L. mandshuricus (smaller species) are inconsistent with this hypothesis. Species abundance and body size, though, might have been different in ancestral populations. The present study indicates that at least some of the hybridization events are recent, therefore, the forces driving the hybridization among Chinese Lepus species may include species abundance asymmetries and mating preferences.
The hybridizations among Chinese Lepus species might be a consequence of range expansion caused by the Pleistocene glaciations. Melo-Ferreira et al.  estimated the timing of mtDNA introgression from L. timidus into L. granatensis in Iberia at 33 000-35 000 years and suggested that it might have been favored by the climatic fluctuation during the Pleistocene glacial phenomena. Our latest study also strongly indicates that one of the Chinese Lepus species (Lepus yarkandensis) has undergone repeated population reductions and expansions during the Pleistocene glacial and interglacial periods, respectively . These population size fluctuations may have contributed to the extensive gene flow among Lepus yarkandensis populations. Similarly, the alternation of glacial and interglacial periods during the Pleistocene likely contributed to the conditions necessary for of the overlapping ranges of different Chinese hare species, resulting in the observed introgressive hybridization. Morgan et al.  and Melo-Ferreira et al.  drew similar inferences for both mosquitoes in Southeast Asia and European hares, respectively. Further support for the effect of glaciation comes from the estimation of some of the introgression events among the Chinese hare lineages as occurring during the Pleistocene (i.e. 0.414-0.988 MYA). Simulation studies have also suggested that range expansions could lead to asymmetric introgression from the resident species toward the invading species [[64–66]].
In addition to the above factors, the widespread and complex patterns of reticulate evolutionary events among the Chinese Lepus lineages may also reflect the generally low levels of interspecific genetic differentiation. Previous studies have suggested that Lepus experienced rapid radiation [67, 68]. The lack of chromosomal structural changes supports such a model of relatively recent, and rapid, diversification in this clade [69, 70]. Low levels of genetic diversification between species could result in relatively fertile F1 hybrids providing a bridge for further hybridization and thus introgression. Alternatively, low levels of genetic diversification between species may indicate that speciation processes of some Chinese Lepus species are not complete and there have been continuous gene flow among them. For example, the admixture of MGF gene between L. sinensis and L. capensis could be the consequence of an incomplete reproductive isolation between the two species. Further definition of the introgressive system characterizing Chinese Lepus species should provide a clearer resolution of the evolutionary factors that have resulted in this reticulate complex.
Taxonomic classification of Chinese hares
Based on the mtDNA and nuclear gene sequence analyses, the present study not only suggests that introgressive hybridization among Chinese Lepus species contributed to the current taxonomic confusion, but also provides insights into the previously obscure taxonomic classification of Chinese hares. Our mtDNA and nuclear DNA analyses supported the specific status of eight hare species, including L. hainanus, L. comus, L. oiostolus, L. yarkandensis, L. timidus, L. mandshuricus, L. capensis and L. sinensis. This finding is consistent with the morphological taxonomy of Pan et al. . The division of L. capensis into two species (L. tibetanus and L. tolai), as suggested by Hoffmann & Smith , was not favored in the present study.
In the mtDNA analysis of L. mandshuricus, nine specimens were found to harbor L. timidus mtDNA haplotypes and three specimens contained L. sinensis mtDNA haplotypes, suggesting the absence of a diagnostic L. mandshuricus mtDNA lineage (Figure 1); however, the nuclear gene (MGF) analysis supported the monophyly of L. mandshuricus, which is consistent with its specific status (Figure 2). Furthermore, Liu et al. , using skull characters and nuclear gene sequences, also supported the specific status of L. mandshuricus, therefore, the discrepancy between the mtDNA and nuclear gene trees is most likely due to mtDNA introgression from L. timidus and L. sinensis into L. mandshuricus. This finding reflects the importance of including nuclear markers for systematic studies of hares [[71–73]].
On the mtDNA tree, L. capensis and L. timidus formed sister taxa and have the lowest Cyt b sequence divergence (2.2%) among the species analyzed in this study (Figure 1, Additional file 1). The Cyt b sequence divergence of 7.6% between L. capensis and L. sinensis (Additional file 1), and the presence of L. capensis group and L. sinensis group on nuclear gene tree (Figure 2) supported the specific status of L. capensis and L. sinensis. As we discussed above, our analyses led to the inference that L. timidus and L. sinensis mitochondrial genomes had introgressed into L. capensis and that frequent nuclear DNA introgression had occurred between L. capensis and L. sinensis. These observations are inconsistent with those of Wu et al. , who concluded that L. capensis either does not exist in China as a unique taxon or has been replaced by L. timidus. In contrast, our data suggest the occurrence of continuous introgression from L. timidus into L. capensis.
In the present study, our mtDNA tree identified a new mtDNA lineage that might represent a previously unidentified hare species from Xinjiang Province. This is evidenced by the fact that the newly identified mtDNA lineage (L. capensis-2) was closely related to L. yarkandensis on the mtDNA tree (Figure 1) and the Cyt b sequence divergence (2.5%) between L. capensis-2 and L. yarkandensis is higher than that (2.2%) found between L. capensis and L. timidus (Additional file 1). This finding is consistent with the inference of Wu et al.  that there might be two or more new hare species in this area of China; however, the MGF haplotypes of the L. capensis-2 samples were divided into two different groups: (i) two (Hap13 and Hap 14) were nested in the L. yarkandensis clade, with one (Hap13) shared with L. yarkandensis; (ii) five (Hap6, 15, 17, 19, 21) fell into the L. capensis group, with two (Hap15 and Hap21) shared with L. capensis (Figure 2). In addition, the individuals of L. capensis-2 and L. capensis have no significant morphological differences, therefore, the alternative hypothesis, that this newly identified mtDNA lineage may reflect ancient mtDNA introgression from L. yarkandensis into L. capensis, seems more likely.