Abstract

The hares and rabbits belonging to the family Leporidae have a nearly worldwide distribution and approximately 72% of the genera have geographically restricted distributions. Despite several attempts using morphological, cytogenetic, and mitochondrial DNA evidence, a robust phylogeny for the Leporidae remains elusive. To provide phylogenetic resolution within this group, a molecular supermatrix was constructed for 27 taxa representing all 11 leporid genera. Five nuclear (SPTBN1, PRKCI, THY, TG, and MGF) and two mitochondrial (cytochrome b and 12S rRNA) gene fragments were analyzed singly and in combination using parsimony, maximum likelihood, and Bayesian inference. The analysis of each gene fragment separately as well as the combined mtDNA data almost invariably failed to provide strong statistical support for intergeneric relationships. In contrast, the combined nuclear DNA topology based on 3601 characters greatly increased phylogenetic resolution among leporid genera, as was evidenced by the number of topologies in the 95% confidence interval and the number of significantly supported nodes. The final molecular supermatrix contained 5483 genetic characters and analysis thereof consistently recovered the same topology across a range of six arbitrarily chosen model specifications. Twelve unique insertion-deletions were scored and all could be mapped to the tree to provide additional support without introducing any homoplasy. Dispersal-vicariance analyses suggest that the most parsimonious solution explaining the current geographic distribution of the group involves an Asian or North American origin for the Leporids followed by at least nine dispersals and five vicariance events. Of these dispersals, at least three intercontinental exchanges occurred between North America and Asia via the Bering Strait and an additional three independent dispersals into Africa could be identified. A relaxed Bayesian molecular clock applied to the seven loci used in this study indicated that most of the intercontinental exchanges occurred between 14 and 9 million years ago and this period is broadly coincidental with the onset of major Antarctic expansions causing land bridges to be exposed.

The order Lagomorpha contains two families, the Ochotonidae (pikas) and Leporidae (rabbits and hares). Although Ochotonidae were more diverse during the Miocene, they are currently characterized by a single extant genus, Ochotona, including approximately 26 species (Angermann et al., 1990; Nowak, 1999; Yu et al., 2000). On the other hand, the rabbits and hares are more species rich, with 11 genera and more than 56 species (Angermann et al., 1990). From a paleontological perspective, it is not clear whether the order originated in Asia or North America and unfortunately a more complete lagomorph fossil record for the early and middle Eocene is lacking (Dawson, 1981). However, it has been postulated that the first expansion of Leporidae occurred in North America during the Miocene (Dawson, 1981), a notion supported by at least two recent fossil discoveries suggesting a North American origin for the family (White, 1991; Voorhies and Timperley, 1997).

Although members of the Leporidae have a nearly worldwide distribution, at least eight of the genera (Brachylagus, Pentalagus, Caprolagus, Bunolagus, Poelagus, Romerolagus, Oryctolagus, and Nesolagus) have geographically restricted distributions (Fig. 1). Apart from Nesolagus—which includes two species (Surridge et al., 1999; Can et al., 2001)—all these genera are monotypic. Of the more widely distributed taxa, Sylvilagus contains more than 16 species restricted to North, Central, and South America (Angermann et al., 1990; Chapman et al., 1992; Frey et al., 1997), whereas the African genus Pronolagus includes at least four species (Angermann et al., 1990; Whiteford, 1995; Matthee and Robinson, 1996). The genus Lepus (hares and jackrabbits) is characterized by approximately 26 species and is the only taxon with an almost cosmopolitan distribution (Flux and Angermann, 1990).

Figure 1

Approximate current distribution of 10 leporid genera: I = Oryctolagus; II = Caprolagus; III = Pentalagus; IV = Nesolagus; V = Poelagus; VI = Bunolagus; VII = Pronolagus; VIII = Romerolagus; IX = Brachylagus; X = Sylvilagus. Arrowheads from circles indicate polytypic species with a wider distribution. The distribution of Lepus is nearly cosmopolitan and not indicated on the map (see text for detail).

Figure 1

Approximate current distribution of 10 leporid genera: I = Oryctolagus; II = Caprolagus; III = Pentalagus; IV = Nesolagus; V = Poelagus; VI = Bunolagus; VII = Pronolagus; VIII = Romerolagus; IX = Brachylagus; X = Sylvilagus. Arrowheads from circles indicate polytypic species with a wider distribution. The distribution of Lepus is nearly cosmopolitan and not indicated on the map (see text for detail).

Resolving leporid evolutionary relationships, particularly among genera, has proven difficult with conventional phylogenetic approaches. In large part this is due to convergence in anatomical features, an absence of chromosomal synapomorphies, and the saturation of mitochondrial DNA sequences. Of the early phylogenetic attempts, those based on premolar tooth patterns were the most definitive (Dawson, 1958; Hibbard, 1963; Dawson, 1981). The most recent of these incorporated 9 of the 11 leporid genera (Poelagus and Bunolagus were not included). The conclusion by Dawson (1981) was that Nesolagus and Brachylagus were basal sister taxa and the remaining leporid genera were subdivided into three equidistant evolutionary lineages (in other words Nesolagus and Brachylagus were morphologically primitive, and an additional three leporid lineages originated more recently and contemporaneously). The first of the derived lineages comprised Pronolagus and Pentalagus as sister taxa, the second grouped Lepus, Oryctolagus, Caprolagus, and Sylvilagus together, while the third lineage was represented by the monotypic Romerolagus diazi (Dawson, 1981). Subsequent to Dawson's (1981) study, Corbet (1983) examined 21 morphological character states for 22 leporid species representative of all the genera and concluded there is considerable homoplasy in leporid morphology. Likewise, the analysis of partial sequences of two mitochondrial genes (cytochrome b and 12S rRNA; Halanych and Robinson, 1999) failed to fully resolve the phylogenetic relationships among most genera, leading the authors to conclude that most genera originated contemporaneously approximately 12 to 16 MYA. Subsequently, Su and Nei (1999) suggested that this date is closer to 20 MYA. Most recently, cross-species chromosome painting (Robinson et al., 2002) showed that the numerous chromosomal rearrangements that characterize the karyotypes of many of the extant leporid genera are autapomorphic, and these probably arose as a result of a combination of low population numbers and founder events. The failure of the chromosomal data to provide markers that track descent from common ancestry, particularly at the intergeneric level, was thought to reflect the group's rapid radiation, thus mirroring the pattern suggested by the mtDNA sequence data.

In instances where conventional phylogenetic reconstruction methods (mainly based on single evolutionary markers) have failed to recover a robust phylogeny (Irwin et al., 1991; Gatesy et al., 1997; Matthee and Robinson, 1999), a supermatrix approach has proved to be an effective way of obtaining better resolved phylogenetic trees (Baker and DeSalle, 1997; Cognato and Vogler, 2001; Murrell et al., 2001; Buckley et al., 2002; Murphy et al., 2002; Gatesy et al., 2002). This applies even in instances where the evolutionary history has been confounded by a rapid radiation (Matthee et al., 2001). To this end, and prompted by the availability of a suite of unique nuclear DNA intron markers that were developed and used with success to generate a molecular supermatrix for the Cetartiodactyla (Matthee et al., 2001) and the family Bovidae (Matthee and Davis, 2001), we followed this approach with the leporids.

Our aims were threefold. First, we hoped that the greater resolving power of nuclear intron sequences would shed new light on the seemingly intractable problem of rapid cladogenesis in the Leporidae. Secondly, we examined the usefulness of separate data partitions and investigated the topological effects of model-selection criteria on these in a data set that is noteworthy for full representation at the generic level. Thirdly, the leporids constitute an ideal group for biogeographic analyses because of their demographic characteristics and wide distribution. Most methods available to tests different biogeographical scenarios are based solely on single branching patterns, which are inadequate given that both geographic areas and taxa are often characterized by multiple histories and reticulate relationships (Morrone and Crisci, 1995). In the present study we use the leporid phylogenetic tree as reference and employed the parsimony based dispersal-vicariance (DiVa) method of Ronquist (1997) to infer past evolutionary biogeographic scenarios that may have caused the distribution patterns observed in extant leporid taxa. Dating of these scenarios will allow extrapolation to other nonrelated taxa and to this end we included a relaxed Bayesian molecular clock using multiple loci to provide posterior probability divergence estimates for nodes. The latter method not only allows for rate differences among taxa but also provides insights into whether evolutionary rates for two genes change in a correlated fashion.

Materials and Methods

Taxonomic Sampling

Representatives of all 11 leporid genera were included. In all polytypic genera (Nesolagus, Pronolagus, Lepus, and Sylvilagus) we examined multiple species resulting in a total of 25 ingroup taxa (Table 1). Their sequences were analyzed in conjunction with two outgroups of increasing evolutionary distance: Ochotona princeps, a representative of the sister family of leporids, the Ochotonidae, and Tamiasciurus hudsonicus, of the order Rodentia, which, together with the Lagomorpha, constitute the superordinal grouping Glires (as recently confirmed by Murphy et al., 2001).

Table 1

Taxonomic designation and common names of species included in the present study. The number of genes sequenced for each taxon are indicated.

Species Common name Number of genes included 
1. Oryctolagus cuniculus European rabbit 
2. Brachylagus idahoensis Pygmy rabbit 
3. Bunolagus monticularis Riverine rabbit 
4. Romerolagus diazi Teporingo, Volcano rabbit 
5. Pentalagus furnessi Amami rabbit 
6. Pronolagus rupestris Smith's red rock rabbit 
7. Pronolagus saundersiae Hewitt's red rock rabbit 
8. Pronolagus randensis Jameson's red rock rabbit 
9. Pronolagus crassicaudatus Natal red rock rabbit 
10. Nesolagus netscheri Sumatran rabbit 
11. Nesolagus timminsi Annamite striped rabbit 
12. Sylvilagus audubonii Desert cottontail 
13. Sylvilagus floridanus Eastern cottontail 
14. Sylvilagus aquaticus Swamp rabbit 
15. Sylvilagus palustris Marsh rabbit 
16. Sylvilagus obscurus Appalachian cottontail 
17. Sylvilagus nuttallii Mountain cottontail 
18. Poelagus marjorita Bunyoro rabbit 
19. Caprolagus hispidus Hispid hare 
20. Lepus americanus Snowshoe hare 
21. Lepus timidus Mountain hare 
22. Lepus townsendii White-tailed jack rabbit 
23. Lepus californicus Black-tailed jack rabbit 
24. Lepus capensis Cape hare 
25. Lepus saxatilis Scrub hare 
26. Ochotona princeps American pika 
27. Tamiasciurus hudsonicus Red squirrel 
Species Common name Number of genes included 
1. Oryctolagus cuniculus European rabbit 
2. Brachylagus idahoensis Pygmy rabbit 
3. Bunolagus monticularis Riverine rabbit 
4. Romerolagus diazi Teporingo, Volcano rabbit 
5. Pentalagus furnessi Amami rabbit 
6. Pronolagus rupestris Smith's red rock rabbit 
7. Pronolagus saundersiae Hewitt's red rock rabbit 
8. Pronolagus randensis Jameson's red rock rabbit 
9. Pronolagus crassicaudatus Natal red rock rabbit 
10. Nesolagus netscheri Sumatran rabbit 
11. Nesolagus timminsi Annamite striped rabbit 
12. Sylvilagus audubonii Desert cottontail 
13. Sylvilagus floridanus Eastern cottontail 
14. Sylvilagus aquaticus Swamp rabbit 
15. Sylvilagus palustris Marsh rabbit 
16. Sylvilagus obscurus Appalachian cottontail 
17. Sylvilagus nuttallii Mountain cottontail 
18. Poelagus marjorita Bunyoro rabbit 
19. Caprolagus hispidus Hispid hare 
20. Lepus americanus Snowshoe hare 
21. Lepus timidus Mountain hare 
22. Lepus townsendii White-tailed jack rabbit 
23. Lepus californicus Black-tailed jack rabbit 
24. Lepus capensis Cape hare 
25. Lepus saxatilis Scrub hare 
26. Ochotona princeps American pika 
27. Tamiasciurus hudsonicus Red squirrel 

Data Collection

With the exception of N. netscheri and Poelagus for which only museum skin samples were available, total genomic DNA was extracted from ethanol preserved tissue or frozen fibroblast cells using a phenol/chlorophorm/isoamyl alcohol procedure (Hillis et al., 1996). In the case of Poelagus, DNA was extracted twice from dried museum skin using the QIAamp DNA purification kit (Qiagen Ltd.). This was done to reduce the risk of contamination by other DNAs. Both extractions resulted in identical sequencing results. DNA from N. netscheri used in this study correspond to those used by Surridge et al. (1999). Attempts were made to amplify five unlinked nuclear DNA fragments (SPTBN1, PRKCI, THY, TG, and MGF) using the A and B primers from Matthee et al. (2001). All primers are situated in the exon regions of the five conserved genes and generated polymerase chain reaction (PCR) fragments mostly containing faster evolving intron sequences. Some of the lagomorph taxa failed to amplify using the published primer sequences and required development of lagomorph-specific primers (Table 2). To obtain amplicons for all individuals (including the outgroup species), lagomorph-specific primers were used by themselves or in combination with the A and B primers mentioned above. In spite of these efforts, amplification of the DNA originating from museum specimens (Poelagus and Nesolagus) proved problematic in some instances (most probably due to the degraded nature of these DNA preparations).

Table 2

Abbreviations and full names of the gene fragments used in the present investigation. PCR annealing temperature and MgCl2 concentrations are given for each primer pair. Leporid specific forward and reverse primers are also shown.

Gene Locus name Annealing temperature PCR MgCl2 PCR Lagomorph forward primer PCR Lagomorph reverse primer 
MGF Stem cell factor 54 2.5 mM 5′-AAATATCAGTCTTGAATCTTAC 5′-TTTTAGATGAATTACAGTGTCC 
TG Thyroglobulin 55 2.5 mM 5′-GCATTGCAGGACAATGAACCCA 5′-CCACTGTTCATACCACTCGAAG 
PRKCI Protein kinase C, iota 57 2.5 mM 5′-AAACAGATCGCATTTATGCAAT 5′-TGTCTGTACCCAGTCAATATC 
SPTBN1 B-spectrin, nonerythrocytic 1 57 2.5 mM 5′-CTCTGCCCAGAAGTTTGCAAC 5′-TGATAGCAGAACTCCATGTGG 
THY Thyrotropin 56 2.5 mM 5′-CATCAACACCACCATCTGTGC 5′-CACTTGCCACACTTACAGCT 
Gene Locus name Annealing temperature PCR MgCl2 PCR Lagomorph forward primer PCR Lagomorph reverse primer 
MGF Stem cell factor 54 2.5 mM 5′-AAATATCAGTCTTGAATCTTAC 5′-TTTTAGATGAATTACAGTGTCC 
TG Thyroglobulin 55 2.5 mM 5′-GCATTGCAGGACAATGAACCCA 5′-CCACTGTTCATACCACTCGAAG 
PRKCI Protein kinase C, iota 57 2.5 mM 5′-AAACAGATCGCATTTATGCAAT 5′-TGTCTGTACCCAGTCAATATC 
SPTBN1 B-spectrin, nonerythrocytic 1 57 2.5 mM 5′-CTCTGCCCAGAAGTTTGCAAC 5′-TGATAGCAGAACTCCATGTGG 
THY Thyrotropin 56 2.5 mM 5′-CATCAACACCACCATCTGTGC 5′-CACTTGCCACACTTACAGCT 

In addition to the nuclear data, all the lagomorph mitochondrial sequences that are available in GenBank were included in our study (Halanych and Robinson, 1997; Gissi et al., 1998; Halanych et al., 1999; Halanych and Robinson, 1999; Surridge et al., 1999; Yamada et al., 2002). To avoid the potential problems associated with missing data and ensure a comprehensive taxon representation (Zwickl and Hillis, 2002; Pollock et al., 2002), we extended the GenBank sequences by adding ∼ 500 additional cytochrome b base pairs for 20 of the 25 species included in our investigation. In the case of the 12S rRNA gene, we provide new sequence data for a further nine species. The standard mammalian primers were used for the mtDNA amplifications and sequencing (Pääbo and Wilson, 1988; Kocher et al., 1989; Irwin et al., 1991; Allard and Honeycutt, 1992).

PCR cycling details followed those described in Matthee et al. (2001). PCR products were purified using the QIAquick PCR purification kit (Quagen Ltd.) and cycle sequenced using BigDye terminator chemistry (Applied Biosystems). Sequencing cocktails were cleaned using Centrisep spin columns (Princeton Separations) and the products were analyzed on a 3100 AB automated sequencer. In most instances, both strands were sequenced to improve the accuracy of base calling and the authenticity of all data generated was determined using GenBank BLASTN searches. A limited number of heterozygous changes occurred in the nuclear fragments and these were designated as one of the two possible states based on their frequency in the other individuals. None of these inferences resulted in the addition of synapomorphic characters and sequences from the 135 nuclear and 29 new mtDNA fragments were deposited in GenBank (THY: AY292663–AY292688; 12S rRNA: AY292689–AY292714; Cyt b: AY292715–AY292738; MGF: AY292739–AY292764; PRKCI:AY292765–AY292791; SPTBN1: AY292808–AY292832; TG: AY292833–AY292859).

Sequence Alignment

Sequences were aligned manually. For the nuclear genes we used the conserved exon sequences to anchor the alignments. Few insertion-deletion events were observed among members of the ingroup and for these, gaps were introduced to maintain the alignment. None of the nuclear insertion-deletions rendered the ingroup alignment problematic; ambiguous alignment was restricted to outgroup comparisons only. Although alternative alignments in these areas are also likely, our investigations showed that minor changes to the outgroup alignment did not significantly alter the relationships among the ingroup taxa (data not presented). A small region of the 12S rRNA sequence (the loop between stems 17 and 18 as defined by Springer et al., 1995) could not be aligned unambiguously and this section was excluded from all analyses. All alignment gaps were treated as missing characters. Insertion-deletion events longer than 2 bp were mapped onto the tree obtained from the analysis of the supermatrix.

Phylogenetic Approach

All fragments were analyzed singly using parsimony (MP) and maximum likelihood (ML) in PAUP v4.0b10 (Swofford, 2000) and Bayesian inference (BI) as implemented in MrBayes v.3.0b3 (Huelsenbeck and Ronquist, 2001). MP analyses were based on heuristic searches with 100 random additions of taxa and tree bisection and reconnection (TBR) branch swapping. Halanych and Robinson (1999) studied the effects of differential weighting on the mtDNA data and showed that the reduction of noise had little effect on improving the leporid mtDNA topologies (also see Källersjö et al., 1999; Cognato and Vogler, 2001). As a result, and given the low prevalence of multiple hits in the nuclear data, no differential weighting of transitions and transversions was applied and no special character weighting scheme was used for each gene partition. One thousand parsimony bootstrap replicates were performed; due to time constraints we restricted the maximum number of trees saved during each replicate to 1000.

ML analyses were performed specifying the optimal model as reflected by the LRT as well as by the AIC criteria implemented in Modeltest v3.06 (Posada and Crandall, 1998). If different optimal models were selected by these two methods, ML topologies were compared for congruence/conflict. We performed 100 ML bootstrap replicates on each topology. For consistency we presented only the ML nodal support using the optimal model selected based on the LRT criteria.

Posterior probabilities for the BI analysis were determined by running one cold and three heated chains for 1,000,000 generations. Data sets containing more than one gene fragment were analyzed in a partitioned manner to allow the selection of different optimal models for each partition (Ronquist and Huelsenbeck, 2003). The optimal model, selected under LRT criteria, was specified as prior for each partition. Each analysis was repeated twice and trees were sampled every 100 generations. Stability was determined using the sump command in MrBayes and the first 20,000 generations were excluded as burn-in each time. Posterior probabilities for nodes were based on the remaining 19,600 topologies.

Data were analyzed by each of the three methods of phylogenetic inference using the following partitions: (a) each gene separately, (b) a combined mtDNA data set, (c) a combined nuclear DNA data set and finally, (d) a supermatrix including all data from all genes. We tested the utility of each partition by comparing the topologies obtained from the partitioned analyses (a–c) to that of the supermatrix (d). The supermatrix topology was chosen on the grounds that 10 nodes were consistently recovered from the supermatrix data set regardless of the method of analysis (the nodes were labeled A to J) and that each identified a monophyletic lineage comprising at least two genera. In addition, most of these nodes had strong bootstrap support (> 70%) and posterior probability values (> 95%).

To investigate possible conflict among the different gene fragments, we used a parsimony approach (ILD as implemented in PAUP v4.0b10; Farris et al., 1995; Baker and De Salle, 1997) and to determine phylogenetic signal, a Bayesian approach (Buckley et al., 2002). For the ILD test, statistical support for incongruence was estimated using 1000 repartitions and a significance value of P < 0.05 was assumed. For the Bayesian approach, the number of unique topologies in the 95% posterior probability was used. In the case of the latter, there is an inverse relationship between the resolving power of the data and the number of trees present in the 95% posterior interval (Buckley et al., 2002). To explore the effects of model selection on the combined maximum likelihood analyses, we repeated the ML analyses of our supermatrix using an arbitrary selection of models of varying parameter richness ranging from zero to 10 (JC = 0; K2P = 1; F81 = 3; HKY85 = 4; GTR = 8; and GTR + I + Γ = 10).

Biogeographic Analyses

For the DiVa analyses, we used the phylogenetic tree based on the supermatrix data set. We defined five geographic unit areas following the current distribution of the leporids and the pika (as the outgroup): North America (A), South America (B), Asia (C), Europe (D), and Africa (E). The current distribution of each genus was coded as present or absent in each of the five unit areas. Although we can also regard Asia and Europe as a single unit area (Eurasia), treating them as such did not affect the main findings because it only affected one node in the topology (the continental dispersal of Oryctolagus/Caprolagus ancestor from Asia to Europe). In our analyses, we set the “maxareas” option of the “optimize” command to two, thereby constraining the ancestral distribution of genera (either North America or Asia) to obtain a narrower estimate of the likely ancestral distributions at each node. This decision was based on the fact that an Asian and/or North American origin for the group is consistent with the fossil record and since there is no evidence to suggest an African origin, we did not take this latter possibility into account (Voorhies and Timperley, 1997).

Molecular Clock

In an attempt to date the evolution of the Leporidae in a biogeographic context, we employed a molecular clock. First, to determine whether rate differences were present among leporid lineages, the relative rates of evolutionary change were compared using RRTree specifying the K2P model of substitution (Robinson-Rechavi and Huchon, 2001). The closest evolutionary outgroup (O. princeps) was used as the reference taxon. This test was chosen because it takes the phylogenetic branching pattern into account and corrects for sampling imbalances (for example multiple species included forLepus, Sylvilagus, Nesolagus, and Pronolagus compared to the remaining monotypic genera). It is important to realize, however, that using an equal rates model (K2P available in the software) when a more complex scenario is present can result in biases favoring a constant molecular clock. This is because long branches will be underestimated disproportionately more than short branches using K2P and this will cause the test statistic (branch-length difference) to be underestimated. Despite using a simple model of evolution, significant rate differences were observed among leporid lineages and it seems reasonable to suggest that the error introduced by the use of an inappropriate model did not affect our result.

Given the absence of a constant molecular clock for the leporids, we employed the recently developed relaxed Bayesian methods for multilocus data (Thorne et al. 1998; Thorne and Kishino, 2002). Four software applications were used: First, in the package PAML v3.14 (Yang, 1997) we used BASEML to create the maximum likelihood output files needed for the MULTIDIVERGENCE package developed by Thorne et al. (1998) and Kishino et al. (2001). We used PAML2MODELINFO to write “model files” needed to estimate the variance-covariance matrix for all seven DNA fragments independently. The latter procedure was performed using ESTBRANCHES. Finally, the output files from ESTBRANCHES were employed in MULTIDIVTIME to estimate the prior and posterior distribution of divergence dates. In addition, these output files from MULTIDIVTIME also incorporates the rank correlation coefficient for each pair of genes and the approximate P-value for evaluating the null hypothesis that the pair of genes change rate independently. The Markov chain Monte Carlo analyses involved an initial burn-in (100,000 cycles), after which the Markov chain was sampled 10,000 times every 100th cycle. As priors we adopted 40 million years (SD = 40 million years) between the tip and the root and 0.003 (SD = 0.003) substitutions per site per million years for the rate at the root node. The value of the rate of the evolution at the root was varied and it was found that quite large changes to the root rate had little influence on the estimates. Changing the constraints on the nodes affected the estimates much more. The Brownian rate parameter was set to the default value in the program.

Tamiasciurus was used as outgroup and to obtain reasonably narrow posterior distributions for divergence times we have incorporated six time constraints from the literature/fossil record. The first pair was the 20 MYA lower and 40 MYA upper divergence between the ochotonid and the leporid lineages (Dawson, 1981; Diersing, 1984; also see Springer et al., 2003); the second was the origin of modern leporids with a lower at 12 MYA (Voorhies and Timperley, 1997) and the upper 20 MYA (Su and Nei, 1999). The third pair involved the 4 MYA lower and 6 MYA upper divergence of the genus Lepus (Yamada et al., 2002). The highest possible divergence time for the ingroup was set at 60 million years. Finally, because the RRTree analyses indicated significant rate heterogeneity among lineages particularly when the mtDNA data were included, we examined the effect of these on the initial analyses. The data were reanalyzed with MULTIDIVTIME excluding all fragments indicating large rate differences. For the latter we kept all priors the same as for the complete data set.

Results

Sequence Data

The combined molecular data set comprised a total of 5483 nucleotides of which 1882 base pairs originated from two genes situated on the mitochondrial genome (Table 3). The remainder of the supermatrix sequences were derived from five nuclear DNA segments and these varied in length from 556 bp in TG to 865 bp in SPTBN1 (excluding some autapomorphic insertions, see below). Except for Poelagus for which only two nuclear genes were sequenced, and Nesolagus where data for cytochrome b and one species' SPTBN1 sequence are lacking, our data matrix is complete (Table 1).

Table 3

Characteristics of the gene segments used in this investigation. Gene names correspond to those in Table 2. The number of taxa used in each analysis is presented and the values for each segment include the total number of characters (including alignment gaps) and the number of parsimony informative and variable characters. The number of equally parsimonious trees, the retention index value (RI), and the number of trees in the 95% posterior probability analyses are given for each data set (see text for detail).

Gene Total number of taxa Total number of characters Parsimony informative characters Variable characters MP tree length RI value Number of equally parsimonious trees Number of trees in 0.95 posterior intervals 
MGF 26 764 101 295 388 0.869 2171 
PRKCI 27 557 76 235 315 0.826 96 4623 
SPTBN1 25 865 116 336 461 0.839 3075 
THY 26 859 82 212 269 0.843 3593 
TG 27 556 90 253 359 0.718 540 4251 
All NuDNA 27 3601 465 1331 1823 0.800 24 
Cyt b 25 1140 429 513 2181 0.440 1023 
12S rRNA 26 742 175 262 729 0.507 4537 
All MtDNA 26 1882 604 775 2929 0.450 2347 
Supermatrix 27 5483 1069 2106 4775 0.560 
Gene Total number of taxa Total number of characters Parsimony informative characters Variable characters MP tree length RI value Number of equally parsimonious trees Number of trees in 0.95 posterior intervals 
MGF 26 764 101 295 388 0.869 2171 
PRKCI 27 557 76 235 315 0.826 96 4623 
SPTBN1 25 865 116 336 461 0.839 3075 
THY 26 859 82 212 269 0.843 3593 
TG 27 556 90 253 359 0.718 540 4251 
All NuDNA 27 3601 465 1331 1823 0.800 24 
Cyt b 25 1140 429 513 2181 0.440 1023 
12S rRNA 26 742 175 262 729 0.507 4537 
All MtDNA 26 1882 604 775 2929 0.450 2347 
Supermatrix 27 5483 1069 2106 4775 0.560 

Based on the number of variable characters available for phylogenetic inference, both the mtDNA and nuclear DNA data sets appeared to be potentially useful (37% of the nuclear DNA and 41% of the mtDNA characters were variable; Table 3). There were, however, marked differences in their utility under parsimony reconstruction when the percentage synapomorphic characters were compared for each (13% synapomorphic characters for the nuclear DNA data and 32% for the mtDNA data). Given the rate of nucleotide substitution, the levels of homoplasy (as assessed by the RI values) were much higher in the mtDNA data (Table 3).

Genetic distances based on the GTR + I + Γ for the combined mtDNA data among species within the same genus ranged from 0.031 substitutions per site between L. townsendii and L. timidus to 0.249 between P. crassicaudatus and P. saundersiae. The combined mtDNA values between genera ranged from 0.121 between N. netscheri and L. townsendii to 0.444 between P. crassicaudatus and Brachylagus. For the combined (more conservative) nuclear fragments, intrageneric values ranged from 0.0024 between L. townsendii and L. timidus to 0.0186 between S.palustris and S. floridanus, whereas intergeneric values were from 0.030 between Pentalagus and Oryctolagus to 0.0691 between P. rupestris and Brachylagus.

Insertions and Deletions

The introns of certain taxa were characterized by unique insertions-deletions longer than 50 bp. For example, all six Lepus had a ∼ 220-bp insert in the intron of THY. Limited interspecific variation was evident in this intron but it was retained in the analyses. A BLAST search of the insertion indicated no significant similarity to any sequence currently available in GenBank. A unique SPTBN1 insertion and THY deletion were detected in Caprolagus (in excess of 320bp in both instances). Interestingly, the SPTBN1 insert showed significant similarity (P < 0.001) to a type C repeat found in the genomes of Oryctolagus,Sylvilagus, and Lepus (not shown). Because the SPTBN1 insert was restricted to a single taxon it was excluded from our phylogenetic analyses.

Twelve shorter phylogenetically informative indels (2 to 18 bp in length) were detected in the data. None of these were homoplasious on the supermatrix tree. Single base pair deletions were not scored because of the potential for sequencing errors while overlapping indels of different lengths were excluded due to possible ambiguities in alignment.

Model Selection

Using the LRT and AIC criteria, the same optimal model of nucleotide evolution was selected for three partitions (12S rRNA, combined mtDNA, all data; Table 4). For the remaining seven comparisons where different models of evolution were selected, the respective ML analyses using alternative models resulted in no conflicting nodes with high bootstrap support. For consistency we present the ML topologies based on LRT criteria only (also see Posada and Crandall, 2001).

Table 4

Comparison between optimal models selected by Modeltest under the LRT and AIC criteria for each gene fragment. Gene names are indicated together with the likelihood score for each fragment.

 LRT AIC 
 
 

 
Gene Model 1-score Model 1-score 
MGF HKY + Γ (1.928) 2941.7345 GTR + Γ (1.850) 2932.4725 
PRKCI HKY + Γ (1.472) 2245.5851 TVM + Γ (1.346) 2232.2594 
SPTBN1 K81uf + Γ (1.135) 3441.9936 GTR + Γ (1.185) 3438.3744 
TG K80 + Γ (0.915) 2512.9290 SYM + Γ (0.909) 2507.6543 
THY HKY + Γ (0.797) 2569.7963 TrN + I (0.4208) 2569.5583 
All nuclear TrN + Γ (1.049) 14108.8320 GTR + Γ (1.075) 14088.16180 
Cyt b GTR + I (0.512) + Γ (1.162) 9949.1490 K81uf + I (0.512) + Γ (1.118) 9950.6119 
12S rRNA TrN + I (0.4742) + Γ (0.5588) 4256.2222 TrN + I (0.4742) + Γ (0.5588) 4256.2222 
All mtDNA GTR + I (0.499) + Γ (0.827) 14335.7950 GTR + I (0.499) + Γ (0.827) 14335.7950 
Supermatrix GTR + I (0.312) + Γ (0.459) 29734.4994 GTR + I (0.312) + Γ (0.459) 29734.4994 
 LRT AIC 
 
 

 
Gene Model 1-score Model 1-score 
MGF HKY + Γ (1.928) 2941.7345 GTR + Γ (1.850) 2932.4725 
PRKCI HKY + Γ (1.472) 2245.5851 TVM + Γ (1.346) 2232.2594 
SPTBN1 K81uf + Γ (1.135) 3441.9936 GTR + Γ (1.185) 3438.3744 
TG K80 + Γ (0.915) 2512.9290 SYM + Γ (0.909) 2507.6543 
THY HKY + Γ (0.797) 2569.7963 TrN + I (0.4208) 2569.5583 
All nuclear TrN + Γ (1.049) 14108.8320 GTR + Γ (1.075) 14088.16180 
Cyt b GTR + I (0.512) + Γ (1.162) 9949.1490 K81uf + I (0.512) + Γ (1.118) 9950.6119 
12S rRNA TrN + I (0.4742) + Γ (0.5588) 4256.2222 TrN + I (0.4742) + Γ (0.5588) 4256.2222 
All mtDNA GTR + I (0.499) + Γ (0.827) 14335.7950 GTR + I (0.499) + Γ (0.827) 14335.7950 
Supermatrix GTR + I (0.312) + Γ (0.459) 29734.4994 GTR + I (0.312) + Γ (0.459) 29734.4994 

ML analyses of our supermatrix employing a variety models (JC, K2P, F81, HKY85, GTR, and GTR + I + Γ), recovered identical branching patterns in each instance. Although different models can sometimes result in significantly different phylogenetic estimates (Kelsey et al., 1999; Whelan et al., 2001; Sullivan and Swofford, 2001), in our case, the arbitrarily chosen models had little effect on the supermatrix estimate of the phylogeny (also see Bollback, 2002; Buckley et al., 2002; Dowton and Austin, 2002).

Independent Nuclear DNA Topologies

Apart from consistent support for the monophyly of the Leporidae (node A) and the monophyly of the clade containing Nesolagus and Pronolagus (node B), the analyses of each nuclear locus separately failed to recover more than 50% of the nodes (Table 5). Similarly, the mtDNA data yielded little phylogenetic structure, and the independent analyses recovered on average between 13% and 33% of the intergeneric nodes (Table 5). Importantly, however, no significant posterior probabilities or bootstrap values were recovered for the nodes where the independent gene topologies were incongruent. This failure to recover intergeneric associations is probably due to the limited number of phylogenetically informative characters present when considering the genes singly. The lack of resolution is corroborated by the high number of trees recovered in the 95% posterior interval (> 2100 each time; Table 3).

Table 5

Congruence between the topologies resulting from the analysis of each gene fragment compared to the topology derived from the supermatrix (Fig. 3). Gene abbreviations correspond to those in Table 2 and the methods of analyses are indicated for each data partition. The nodes (A–J) correspond to those shown on Figure 3 and reflect the resolution among Leporid genera. “X” indicates that the node was not recovered and “??” indicate missing data for Poelagus and Nesolagus (see text for detail). “√” indicates the presence of a node and the values above each represent the statistical support. The percentage of times a node was obtained, the percentage congruence among methods, and the percentage of time the nodes were statistically supported (> 70% bootstrap; > 95% posterior probabilities) are also indicated.

  Node  
  
 
 
Gene Analyses % Congruence 
MGF MP √100 √86 √82 √56 ?? √85 √86 67 
 ML √83 √83 √85 √62 ?? √91 √81 67 
 BI √98 √99 √99 √89 ?? √100 √100 67 
PRKCI MP √100 √55 √85 √96 √54 50 
 ML √67 √67 √68 √ < 50 √88 √ < 50 60 
 BI √95 √86 √95 √99 √95 50 
SPTBN1 MP √100 √92 √71 √ < 50 ?? √71 √56 67 
 ML √86 √100 ?? √70 √55 44 
 BI √99 √100 ?? √84 √91 44 
THY MP √100 √79 ?? √ < 50 √54 √63 56 
 ML √96 √83 ?? √63 √52 √69 56 
 BI √100 √80 ?? √65 33 
TG MP √64 √86 √ < 50 √ < 50 40 
 ML √64 √70 √ < 50 30 
 BI √59 √93 √81 30 
 %Obtained 100 100 47 40 50 67 40 47 0 13  
 %Supported 73 67 40 0 50 40 0 20 0 0  
Nuc total MP √100 √91 √84 √66 √90 √99 √96 √98 80 
 ML √100 √96 √86 √57 √87 √99 √94 √98 80 
 BI √100 √100 √100 √89 √100 √100 √100 √100 80 
 %Obtained 100 100 100 100 100 100 100 100 0 0  
 %Supported 100 100 100 0 100 100 100 100 0 0  
12S rRNA MP √99 √ < 50 ?? 22 
 ML √99 √55 ?? 22 
 BI √100 √73 ?? √88 33 
Cyt b MP √93 ?? ?? 13 
 ML √68 ?? ?? 13 
 BI √95 ?? ?? 13 
 %Obtained 100 100 0 0 ?? 0 17 0 0 0  
 %Supported 83 0 0 0 ?? 0 0 0 0 0  
Mt total MP √100 ?? 11 
 ML √100 √ < 50 ?? √ < 50 √ < 50 44 
 BI √100 √51 ?? √98 √55 44 
 %Obtained 100 66 0 0 ?? 0 33 0 33 66  
 %Supported 100 0 0 0 ?? 0 33 0 0 0  
  Node  
  
 
 
Gene Analyses % Congruence 
MGF MP √100 √86 √82 √56 ?? √85 √86 67 
 ML √83 √83 √85 √62 ?? √91 √81 67 
 BI √98 √99 √99 √89 ?? √100 √100 67 
PRKCI MP √100 √55 √85 √96 √54 50 
 ML √67 √67 √68 √ < 50 √88 √ < 50 60 
 BI √95 √86 √95 √99 √95 50 
SPTBN1 MP √100 √92 √71 √ < 50 ?? √71 √56 67 
 ML √86 √100 ?? √70 √55 44 
 BI √99 √100 ?? √84 √91 44 
THY MP √100 √79 ?? √ < 50 √54 √63 56 
 ML √96 √83 ?? √63 √52 √69 56 
 BI √100 √80 ?? √65 33 
TG MP √64 √86 √ < 50 √ < 50 40 
 ML √64 √70 √ < 50 30 
 BI √59 √93 √81 30 
 %Obtained 100 100 47 40 50 67 40 47 0 13  
 %Supported 73 67 40 0 50 40 0 20 0 0  
Nuc total MP √100 √91 √84 √66 √90 √99 √96 √98 80 
 ML √100 √96 √86 √57 √87 √99 √94 √98 80 
 BI √100 √100 √100 √89 √100 √100 √100 √100 80 
 %Obtained 100 100 100 100 100 100 100 100 0 0  
 %Supported 100 100 100 0 100 100 100 100 0 0  
12S rRNA MP √99 √ < 50 ?? 22 
 ML √99 √55 ?? 22 
 BI √100 √73 ?? √88 33 
Cyt b MP √93 ?? ?? 13 
 ML √68 ?? ?? 13 
 BI √95 ?? ?? 13 
 %Obtained 100 100 0 0 ?? 0 17 0 0 0  
 %Supported 83 0 0 0 ?? 0 0 0 0 0  
Mt total MP √100 ?? 11 
 ML √100 √ < 50 ?? √ < 50 √ < 50 44 
 BI √100 √51 ?? √98 √55 44 
 %Obtained 100 66 0 0 ?? 0 33 0 33 66  
 %Supported 100 0 0 0 ?? 0 33 0 0 0  

Combining Data

The ILD analyses revealed that in only two pairwise comparisons (involving SPTBN1 versus THY and MGF versus 12S rRNA) did the two respective genes exhibit greater intergenic incongruence than two random partitions of a homogeneous data set (Table 6). In the former this involved the placement of Lepus and Romerolagus where THY indicated a sister taxon relationship between the two genera and the Pronolagus/Nesolagus clade, whereas the SPTBN1 topology suggests that Romerolagus and Lepus are basal members of a clade containing Brachylagus, Sylvilagus, Oryctolagus, Bunolagus, Pentalagus, and Caprolagus. In the latter case, MGF resulted in a tree nearly identical to the supermatrix topology, whereas the 12S rRNA data resulted in a topology that, amongst others, clustered Oryctolagus with the Nesolagus/Pronolagus clade, Brachylagus as the most basal leporid, and Caprolagus, Romerolagus, and Pentalagus as sister taxa to Sylvilagus.

Table 6

Incongruent length differences (ILDs) and rank correlation coefficients between evolutionary rates of genes for the pairwise comparisons among the seven fragments included in this study. The values in normal text above the diagonal represent the number of extra steps introduced by combining the partitions and the significance values (P) are given below. The values in italics above the diagonal represent the rank correlation coefficient and the italic values below the diagonal represent the proportion of times that the observed rank correlation coefficient was equaled or exceeded when the distribution for the correlation statistic was approximated under the null hypothesis of independent rate changes among the two genes.

Genes SPTBN1 THY MGF PRKC1 TG 125 rRNA Cyt b 
SPTBN1 — 12 13 
  0.07 0.63 0.43 0.49 0.24 −0.24 
THY 0.04* — 13 15 
 0.91  −0.05 0.53 0.06 −0.35 −0.29 
MGF 0.07 0.09 — 19 14 
 0.14 0.98  0.37 0.28 0.65 −0.04 
PRKCI 0.17 0.21 0.59 — 12 
 0.47 0.35 0.59  0.35 −0.09 −0.51 
TG 0.14 0.06 0.18 0.62 — 13 11 
 0.25 0.86 0.53 0.47  0.13 0.27 
12S rRNA 0.23 0.32 0.01* 0.44 0.24 — 20 
 0.48 0.99 0.030.81 0.80  0.36 
Cyt b 0.71 0.74 0.67 1.0 0.94 0.57 — 
 0.89 0.95 0.71 0.99 0.91 0.21  
Genes SPTBN1 THY MGF PRKC1 TG 125 rRNA Cyt b 
SPTBN1 — 12 13 
  0.07 0.63 0.43 0.49 0.24 −0.24 
THY 0.04* — 13 15 
 0.91  −0.05 0.53 0.06 −0.35 −0.29 
MGF 0.07 0.09 — 19 14 
 0.14 0.98  0.37 0.28 0.65 −0.04 
PRKCI 0.17 0.21 0.59 — 12 
 0.47 0.35 0.59  0.35 −0.09 −0.51 
TG 0.14 0.06 0.18 0.62 — 13 11 
 0.25 0.86 0.53 0.47  0.13 0.27 
12S rRNA 0.23 0.32 0.01* 0.44 0.24 — 20 
 0.48 0.99 0.030.81 0.80  0.36 
Cyt b 0.71 0.74 0.67 1.0 0.94 0.57 — 
 0.89 0.95 0.71 0.99 0.91 0.21  
*

Statistical significance at P < 0.05.

Combining the nuclear data dramatically increased the resolving power of the analyses (Table 3; Fig. 2a). For example, 80% of the nodes defining generic associations were retrieved (mostly with significant statistical support); in contrast, 30% to 67% of the nodes were recovered using individual gene fragments (Table 5). Similarly, in the Bayesian analyses the number of trees in the 95% posterior probability drastically decreased when the data were combined (an average of 3543 trees were recovered in the individual analyses versus four for the combined data set; Table 3). In sharp contrast to the nuclear genes, combining the two mtDNA fragments did not result in an increase in phylogenetic resolution either in terms of the number of nodes recovered (Table 5) or the number of trees in the 95% posterior interval (Table 3). In fact, most nodes in the combined mtDNA topology collapsed to an unresolved polytomy (Fig. 2b). We argue that because there were few statistically supported nodes in the mtDNA topologies that were consistent across methods (also see Halanych and Robinson, 1997), the addition of the mtDNA data to the nuclear genes did not significantly alter any of the previously supported nodes derived from the nuclear DNA alone. Strikingly, when the combined nuclear DNA data were considered and a single representative species of each genus was included (all monotypic genera and S. nuttallii, P. randensis, L. saxatilis, andN. netscheri), a single tree was found in the 0.95 posterior interval. The intergeneric branching pattern is fully congruent with the ML topology presented in Figure 3. It therefore seems probable that the larger number of unique trees found with a 95% probability when all taxa are considered is partly due to branch swapping among species within polytypic genera.

Figure 2

Bootstrap concensus topologies of (a) the combined nuclear phylogeny derived from five DNA fragments (SPTBN1, PRKCI, MGF, THY, TG) and (b) a combined mtDNA phylogeny derived from two mtDNA genes (cyt b and 12S rRNA). Values above and below nodes represent the nodal support (MP = parsimony bootstrap; ML = maximum likelihood bootstrap; and BI = Bayesian inference posterior probability) for associations among the 25 ingroup taxa. Values in circles indicate the node numbers corresponding to those in Table 7.

Figure 2

Bootstrap concensus topologies of (a) the combined nuclear phylogeny derived from five DNA fragments (SPTBN1, PRKCI, MGF, THY, TG) and (b) a combined mtDNA phylogeny derived from two mtDNA genes (cyt b and 12S rRNA). Values above and below nodes represent the nodal support (MP = parsimony bootstrap; ML = maximum likelihood bootstrap; and BI = Bayesian inference posterior probability) for associations among the 25 ingroup taxa. Values in circles indicate the node numbers corresponding to those in Table 7.

Figure 3

Supermatrix ML phylogeny derived from all seven gene fragments used in the present study. Arrows indicate unique insertions/deletion supporting associations. Nodes are labeled A to G and correspond to those in Table 5. Values above and below nodes represent the nodal support (MP = maximum parsimony bootstrap; ML = maximum likelihood bootstrap; and BI = bayesian inference posterior probability) for associations among the 25 ingroup taxa. An asterisk indicates that the node was not supported by > 50% bootstrap support.

Figure 3

Supermatrix ML phylogeny derived from all seven gene fragments used in the present study. Arrows indicate unique insertions/deletion supporting associations. Nodes are labeled A to G and correspond to those in Table 5. Values above and below nodes represent the nodal support (MP = maximum parsimony bootstrap; ML = maximum likelihood bootstrap; and BI = bayesian inference posterior probability) for associations among the 25 ingroup taxa. An asterisk indicates that the node was not supported by > 50% bootstrap support.

Supermatrix Topology and the Evolution of the Leporids

Both ML and MP recovered the same topology and bootstrap support for the 10 nodes defining associations among genera ranged from 74% to 92% (Fig. 3). The BI resulted in an identical topology and the posterior probabilities for seven of the nine nodes were significant (Fig. 3). Two major clades were consistently retrieved: one containing Nesolagus, Poelagus, andPronolagus and the second, the remaining genera (Romerolagus, Lepus, Sylvilagus, Brachylagus, Pentalagus, Bunolagus, Caprolagus, and Oryctolagus). In addition to the high bootstrap values (> 79%), three unique nuclear DNA indels also support this early divergence.

Within the first clade, Nesolagus was consistently placed basal to the African sister taxa Poelagus and Pronolagus with some morphological support for these associations. Nesolagus shows a number of primitive features and Corbet (1983) was unable to identify any derived characters that link it to any of the other leporids. Likewise, Ellerman and Morrison-Scott (1951) regarded Poelagus and Pronolagus as congeneric based on cranial characters.

Within the second clade Romerolagus (regarded as being the most morphologically primitive of the leporids; Corbet, 1983) was always placed basal, followed by Lepus. The remainder of the leporid genera formed a clade within which there was good support for the sister taxon relationship between Sylvilagus and Brachylagus (also see Corbet, 1983; Chapman and Ceballos, 1990; Halanych and Robinson, 1997 for similar conclusions). Furthermore, this was supported by a synapomorphic PRKC1 indel of 3 bp length. In addition, there was unambiguous support for the monophyly of a derived group comprising Pentalagus, Bunolagus, Caprolagus, andOryctolagus (Node H; Table 5).

The monophyly of the polytypic genera was a constant feature of the analysis of the supermatrix and these were further supported by at least one synapomorphic indel in each instance (Fig. 3). Apart from Lepus, with numerous unresolved nodes (also see Halanych et al., 1999), the phylogenetic associations among congeneric species were also resolved (Fig. 3). Pronolagus crassicaudatus was always basal followed by P. randensis as the second most basal taxon with P. rupestris and P. saundersiae as the most apical sister taxa (also see Whiteford, 1995). Apart from a single node that is characterised by conflict between the mitochondrial and nuclear data (Fig. 2a, b), the associations among the Sylvilagus species were also congruent with previous mtDNA and some cytogenetic evidence (Robinson et al., 1984; Halanych and Robinson, 1997).

Biogeography

Restricting the maximum number of areas to two for potential ancestral distributions resulted in several possible solutions. The ancestor to the modern leporids had three possible area solutions (either Asia, North America and Asia, or North America and Africa; Fig. 4). The most parsimonious solution, postulating that the ancestor to the leporids was widespread throughout North America and Asia, required a total of nine dispersal events and five vicariance events to account for current distributions (Fig. 4).

Figure 4

Supermatrix topology reflecting the results of optimizations for the ancestral distributions and the inferred geographic pathways from DiVa analyses. All possible ancestral distributions are indicated above the branches at each node, with the most likely explanation in normal text bold. The unit areas correspond to North America = A, South America = B, Asia = C, Europe = D, and Africa = E. The dispersal events are indicated by dashed lines with cross bars indicating the number of events. Vicariance events are indicated by an asterisk. The approximate age (above the branches) of each node is indicated in millions of years and these were obtained from Table 7.

Figure 4

Supermatrix topology reflecting the results of optimizations for the ancestral distributions and the inferred geographic pathways from DiVa analyses. All possible ancestral distributions are indicated above the branches at each node, with the most likely explanation in normal text bold. The unit areas correspond to North America = A, South America = B, Asia = C, Europe = D, and Africa = E. The dispersal events are indicated by dashed lines with cross bars indicating the number of events. Vicariance events are indicated by an asterisk. The approximate age (above the branches) of each node is indicated in millions of years and these were obtained from Table 7.

Rate Heterogeneity Among Lineages and the Molecular Clock

The relative rate test showed no rate differences among lineages in the case of PRKCI and only taxon specific rate differences for three nuclear genes MGF, SPTBN1, and THY. Specifically, MGF showed significant rate difference between S. aquaticus versus S.audubonii, SPTBN1 between S. floridanus and most other lagomorph species, and THY between Oryctolagus versus S. aquaticus and Oryctolagus versus S. palustris. In contrast, TG and both mtDNA genes showed significant rate heterogeneity that involved the majority of species (> 80% of the comparisons).

The relaxed Bayesian analyses indicated that the posterior molecular divergence dates derived from all seven gene fragments had a narrower credibility interval than the priors (Table 7), indicating that most of the molecular dating comes from the data and not the priors (Hassanin and Douzery, 2003). The second round of analyses, excluding the fragments characterized by high levels of rate heterogeneity (12S rRNA, cyt b, and TG), resulted in nearly identical divergence estimates suggesting that, as expected, rate heterogeneity among lineages has little effect on the Bayesian estimation of the divergence times. Interestingly, however, the CredI95% was narrower when the genes showing considerable rate heterogeneity were excluded. This was particularly evident towards the tips of the tree (nodes 1 to 16 Fig. 2a; Table 7). Another prominent result to emerge was the difference in the divergence estimates of node 5 reflecting the clustering of Poelagus with Pronolagus. When all genes were compared to each other, the divergence estimate was 7.98 MYA but when the genes showing rate heterogeneity were excluded (resulting in a single gene being analysed for Poelagus), the divergence estimate decreased to 4.35 MYA. We attribute this large difference to the limited data included for Poelagus (in the four-gene example a single gene PRKCI was included for Poelagus) and the associated problems with the relaxed Bayesian methods when single genes are analysed (Yang and Yoder, 2003). Clearly, the inclusion of multiple genes/fragments will provide more robust estimates of divergence times (Thorne and Kishino, 2002; Yang and Yoder, 2003).

Table 7

Prior and posterior divergence estimates for the nodes defining the associations among all species included in the present study. Node numbers correspond to those indicated on Figure 2a. The time is given in million of years, the standard deviation and 95% credibility interval is given for each node. The first set of values represent the prior and posterior estimates of four nuclear genes where limited rate heterogeneity was evident among taxa while the second set represent the analyses based on all genes, including those showing significant rate heterogeneity among ∼ 80% of the taxa (see text for detail).

 Ages (MYA) 
 
 
 Prior distribution Posterior distribution 
 
 

 
Node Date S.D. 95% CI Date S.D. 95% CI 
SPTBN/THY/PRKCI/MGF 
  1 8.07 (± 5.30) (0.31–19.02) 1.24 (± 0.76) (0.11–3.06) 
  2 3.33 (± 2.88) (0.10–10.57) 1.65 (± 0.69) (0.53–3.21) 
  3 6.57 (± 3.64) (0.97–14.66) 2.30 (± 0.77) (1.08–4.09) 
  4 9.76 (± 3.95) (2.85–17.92) 2.77 (± 0.86) (1.43–4.80) 
  5 13.05 (± 4.06) (5.31–20.82) 4.35 (± 1.54) (2.08–8.05) 
  6 16.23 (± 3.75) (8.57–22.89) 11.42 (± 2.05) (8.10–15.98) 
  7 1.30 (± 1.02) (0.04–3.74) 1.44 (± 0.59) (0.40–2.72) 
  8 2.60 (± 1.19) (0.48–4.93) 1.92 (± 0.63) (0.84–3.32) 
  9 1.98 (± 1.28) (0.08–4.66) 0.73 (± 0.46) (0.05–1.81) 
10 3.91 (± 1.07) (1.48–5.63) 2.48 (± 0.69) (1.29–4.04) 
11 5.16 (± 0.56) (4.09–5.97) 4.79 (± 0.54) (4.03–5.90) 
12 2.60 (± 2.09) (0.10–7.82) 2.25 (± 0.78) (0.97–4.01) 
13 1.76 (± 1.60) (0.05–5.93) 1.24 (± 0.56) (0.32–2.54) 
14 1.79 (± 1.67) (0.05–6.33) 1.95 (± 0.73) (0.72–3.62) 
15 3.53 (± 2.19) (0.48–8.69) 3.66 (± 1.00) (2.01–5.95) 
16 5.31 (± 2.67) (1.34–11.62) 4.21 (± 1.08) (2.43–6.65) 
17 7.12 (± 3.02) (2.42–13.96) 9.03 (± 1.58) (6.43–12.65) 
18 2.35 (± 2.19) (0.06–8.08) 6.87 (± 1.42) (4.55–10.08) 
19 4.65 (± 2.89) (0.71–11.55) 7.82 (± 1.49) (5.35–11.15) 
20 6.88 (± 3.22) (1.95–14.38) 8.64 (± 1.53) (6.11–12.14) 
21 9.09 (± 3.46) (3.72–16.94) 9.95 (± 1.64) (7.31–13.66) 
22 10.99 (± 3.74) (5.28–18.94) 12.80 (± 1.85) (9.88–16.96) 
23 15.20 (± 3.90) (7.72–22.44) 13.61 (± 1.95) (10.60–18.05) 
24 19.32 (± 3.10) (12.92–24.55) 15.28 (± 2.18) (12.20–20.23) 
25 23.89 (± 3.39) (20.12–33.11) 31.73 (± 4.34) (23.31–39.26) 
SPTBN/THY/PRKCI/MGF/TG/Cyt b/12S rRNA 
  1 7.34 (± 4.89) (0.29–17.59) 1.80 (± 0.72) (0.59–3.40) 
  2 3.08 (± 2.73) (0.08–10.10) 2.74 (± 0.63) (1.68–4.13) 
  3 6.17 (± 3.51) (0.91–14.25) 3.93 (± 0.79) (2.58–5.70) 
  4 9.37 (± 3.89) (2.52–17.68) 4.88 (± 0.92) (3.32–6.97) 
  5 12.60 (± 3.95) (5.05–20.21) 7.98 (± 1.75) (4.89–11.67) 
  6 15.74 (± 3.77) (8.04–22.59) 11.30 (± 1.53) (8.73–14.84) 
  7 1.28 (± 1.00) (0.04–3.68) 2.93 (± 0.54) (1.94–4.04) 
  8 2.58 (± 1.17) (0.50–4.89) 3.29 (± 0.56) (2.25–4.42) 
  9 1.89 (± 1.26) (0.08–4.55) 1.22 (± 0.28) (0.74–1.85) 
10 3.87 (± 1.08) (1.51–5.63) 3.70 (± 0.56) (2.66–4.85) 
11 5.17 (± 0.55) (4.09–5.96) 5.16 (± 0.51) (4.14–5.95) 
12 2.96 (± 2.45) (0.11–9.21) 3.38 (± 0.63) (2.30–4.76) 
13 2.04 (± 2.00) (0.05–7.55) 2.18 (± 0.47) (1.38–3.21) 
14 2.05 (± 2.01) (0.06–7.60) 3.22 (± 0.62) (2.17–4.55) 
15 3.96 (± 2.66) (0.50–10.71) 4.52 (± 0.76) (3.20–6.16) 
16 5.85 (± 3.03) (1.34–12.76) 5.33 (± 0.82) (3.91–7.07) 
17 7.69 (± 3.41) (2.52–15.38) 9.22 (± 1.16) (7.23–11.76) 
18 2.31 (± 2.14) (0.06–7.91) 8.47 (± 1.12) (6.52–10.90) 
19 4.80 (± 3.04) (0.76–12.37) 8.94 (± 1.14) (6.98–11.42) 
20 7.09 (± 3.45) (1.89–15.18) 9.44 (± 1.15) (7.46–11.96) 
21 9.40 (± 3.62) (3.93–17.74) 10.28 (± 1.18) (8.29–12.85) 
22 11.16 (± 3.88) (5.42–19.80) 11.80 (± 1.24) (9.79–14.54) 
23 15.14 (± 3.90) (7.60–22.55) 12.77 (± 1.34) (10.73–15.83) 
24 19.25 (± 3.11) (12.93–24.56) 14.00 (± 1.47) (12.08–17.48) 
25 24.02 (± 3.56) (20.13–33.94) 28.96 (± 3.80) (22.44–37.15) 
 Ages (MYA) 
 
 
 Prior distribution Posterior distribution 
 
 

 
Node Date S.D. 95% CI Date S.D. 95% CI 
SPTBN/THY/PRKCI/MGF 
  1 8.07 (± 5.30) (0.31–19.02) 1.24 (± 0.76) (0.11–3.06) 
  2 3.33 (± 2.88) (0.10–10.57) 1.65 (± 0.69) (0.53–3.21) 
  3 6.57 (± 3.64) (0.97–14.66) 2.30 (± 0.77) (1.08–4.09) 
  4 9.76 (± 3.95) (2.85–17.92) 2.77 (± 0.86) (1.43–4.80) 
  5 13.05 (± 4.06) (5.31–20.82) 4.35 (± 1.54) (2.08–8.05) 
  6 16.23 (± 3.75) (8.57–22.89) 11.42 (± 2.05) (8.10–15.98) 
  7 1.30 (± 1.02) (0.04–3.74) 1.44 (± 0.59) (0.40–2.72) 
  8 2.60 (± 1.19) (0.48–4.93) 1.92 (± 0.63) (0.84–3.32) 
  9 1.98 (± 1.28) (0.08–4.66) 0.73 (± 0.46) (0.05–1.81) 
10 3.91 (± 1.07) (1.48–5.63) 2.48 (± 0.69) (1.29–4.04) 
11 5.16 (± 0.56) (4.09–5.97) 4.79 (± 0.54) (4.03–5.90) 
12 2.60 (± 2.09) (0.10–7.82) 2.25 (± 0.78) (0.97–4.01) 
13 1.76 (± 1.60) (0.05–5.93) 1.24 (± 0.56) (0.32–2.54) 
14 1.79 (± 1.67) (0.05–6.33) 1.95 (± 0.73) (0.72–3.62) 
15 3.53 (± 2.19) (0.48–8.69) 3.66 (± 1.00) (2.01–5.95) 
16 5.31 (± 2.67) (1.34–11.62) 4.21 (± 1.08) (2.43–6.65) 
17 7.12 (± 3.02) (2.42–13.96) 9.03 (± 1.58) (6.43–12.65) 
18 2.35 (± 2.19) (0.06–8.08) 6.87 (± 1.42) (4.55–10.08) 
19 4.65 (± 2.89) (0.71–11.55) 7.82 (± 1.49) (5.35–11.15) 
20 6.88 (± 3.22) (1.95–14.38) 8.64 (± 1.53) (6.11–12.14) 
21 9.09 (± 3.46) (3.72–16.94) 9.95 (± 1.64) (7.31–13.66) 
22 10.99 (± 3.74) (5.28–18.94) 12.80 (± 1.85) (9.88–16.96) 
23 15.20 (± 3.90) (7.72–22.44) 13.61 (± 1.95) (10.60–18.05) 
24 19.32 (± 3.10) (12.92–24.55) 15.28 (± 2.18) (12.20–20.23) 
25 23.89 (± 3.39) (20.12–33.11) 31.73 (± 4.34) (23.31–39.26) 
SPTBN/THY/PRKCI/MGF/TG/Cyt b/12S rRNA 
  1 7.34 (± 4.89) (0.29–17.59) 1.80 (± 0.72) (0.59–3.40) 
  2 3.08 (± 2.73) (0.08–10.10) 2.74 (± 0.63) (1.68–4.13) 
  3 6.17 (± 3.51) (0.91–14.25) 3.93 (± 0.79) (2.58–5.70) 
  4 9.37 (± 3.89) (2.52–17.68) 4.88 (± 0.92) (3.32–6.97) 
  5 12.60 (± 3.95) (5.05–20.21) 7.98 (± 1.75) (4.89–11.67) 
  6 15.74 (± 3.77) (8.04–22.59) 11.30 (± 1.53) (8.73–14.84) 
  7 1.28 (± 1.00) (0.04–3.68) 2.93 (± 0.54) (1.94–4.04) 
  8 2.58 (± 1.17) (0.50–4.89) 3.29 (± 0.56) (2.25–4.42) 
  9 1.89 (± 1.26) (0.08–4.55) 1.22 (± 0.28) (0.74–1.85) 
10 3.87 (± 1.08) (1.51–5.63) 3.70 (± 0.56) (2.66–4.85) 
11 5.17 (± 0.55) (4.09–5.96) 5.16 (± 0.51) (4.14–5.95) 
12 2.96 (± 2.45) (0.11–9.21) 3.38 (± 0.63) (2.30–4.76) 
13 2.04 (± 2.00) (0.05–7.55) 2.18 (± 0.47) (1.38–3.21) 
14 2.05 (± 2.01) (0.06–7.60) 3.22 (± 0.62) (2.17–4.55) 
15 3.96 (± 2.66) (0.50–10.71) 4.52 (± 0.76) (3.20–6.16) 
16 5.85 (± 3.03) (1.34–12.76) 5.33 (± 0.82) (3.91–7.07) 
17 7.69 (± 3.41) (2.52–15.38) 9.22 (± 1.16) (7.23–11.76) 
18 2.31 (± 2.14) (0.06–7.91) 8.47 (± 1.12) (6.52–10.90) 
19 4.80 (± 3.04) (0.76–12.37) 8.94 (± 1.14) (6.98–11.42) 
20 7.09 (± 3.45) (1.89–15.18) 9.44 (± 1.15) (7.46–11.96) 
21 9.40 (± 3.62) (3.93–17.74) 10.28 (± 1.18) (8.29–12.85) 
22 11.16 (± 3.88) (5.42–19.80) 11.80 (± 1.24) (9.79–14.54) 
23 15.14 (± 3.90) (7.60–22.55) 12.77 (± 1.34) (10.73–15.83) 
24 19.25 (± 3.11) (12.93–24.56) 14.00 (± 1.47) (12.08–17.48) 
25 24.02 (± 3.56) (20.13–33.94) 28.96 (± 3.80) (22.44–37.15) 

There was generally a negative correlation among the rates of change of genes when the mitochondrial and nuclear fragments were compared; there was a positive correlation when nuclear fragments were compared to each other and when the two mtDNA genes were compared to each other. In one instance the rate differences was significant (MGF versus 12S rRNA; Table 6), but it is also important to realize that the test statistic tends to be positive even when genes evolve independently (Thorne and Kishino, 2002).

Discussion

Resolving Power of the mtDNA Sequence Data

Even with substantially better generic representation and increasing the length of the mitochondrial nucleotide sequences included by Halanych and Robinson (1999) and more recently by Yamada et al. (2002), our investigation did not result in an improvement in the phylogenetic resolution. Saturation, as suggested by low RI values (< 0.5; Table 3), is probably the strongest contributing factor to the weak phylogenetic resolution (Halanych and Robinson, 1999). Moreover, even after applying parameter-rich models that expect multiple substitutions (Whelan et al., 2001), and combining the data from both cytochrome b and 12S rRNA, both ML and Bayesian analyses failed to recover any meaningful phylogenetic structure, or in the case of the latter, to result in a decrease in the number of trees in the 95% credibility interval. In fact, more trees were obtained in the 95% credibility interval for the combined mitochondrial DNA data set than when the cytochrome b gene was analyzed singly. Importantly, however, within the polytypic Pronolagus, Sylvilagus, and Lepus, the mtDNA data were effective in resolving most of the evolutionary relationships among species (also see Halanych and Robinson, 1997; Halanych et al., 1999; Yamada et al., 2002).

Resolving Power of the Nuclear DNA Sequence Data

To date, all single evolutionary marker systems failed to recover any meaningful phylogenetic resolution among the leporid genera (Dawson, 1981; Corbet, 1983; Robinson et al., 1984, 2002; Halanych and Robinson, 1999). This was mainly ascribed to a rapid radiation (Halanych and Robinson, 1997, Robinson et al., 2002), coupled with saturation of the mtDNA data (Halanych, 1997), absence of chromosomal synapomorphies (Robinson et al., 1984, 2002), an incomplete fossil record rendering the interpretation of morphological character evolution problematic (Dawson, 1981), and convergent morphological evolution (Corbet, 1983). The limitation of single markers is further highlighted by our analyses of single nuclear gene fragments—for example when SPTBN1, THY, and TG were analyzed independently, a maximum of 2 out of 10 nodes were consistently recovered with strong support (Table 5). The remaining two nuclear genes performed only marginally better; MGF recovered 5 of 10 intergeneric nodes whereas PRKCI recovered 4 out of 10 intergeneric nodes (Table 5). Moreover, a large number of equally parsimonious trees were found in each instance (Table 3).

In contrast, on combining the five nuclear DNA fragments, 9 of 10 nodes recovered by the supermatrix (Fig. 3) were consistently supported by high bootstrap and posterior probability values (Table 5; Fig. 2a), underscoring the improved resolution often obtained by combining nuclear data (e.g., Baker and DeSalle, 1997; Cognato and Vogler, 2001; Matthee and Davis, 2001; Murrell et al., 2001; Murphy et al., 2001; Buckley et al., 2002). In the case of our data, there was also a considerable decrease in the number of trees obtained in the 95% interval. When the data fragments were analyzed individually, between 2171 and 4623 topologies were found with 95% probabilities (Table 3). On combining the nuclear data only 24 trees were found in the 95% interval.

Resolving Power of the Supermatrix

Combining all seven gene fragments into a single molecular supermatrix resulted in marginally better resolution than that obtained by the combined nuclear matrix (Table 3). Despite the high level of homoplasy present in the mtDNA data, the inclusion of the mitochondrial sequences resulted in improved bootstrap and posterior probability support for one additional intergeneric node (“D,” previously found in three of the five nuclear genes but not supported by either bootstrap or posterior probabilities; Table 5). The sister taxa relationships supported by nodes I and J in the supermatrix tree (O. cuniculus, C. hispidus, B. monticularis, and P. furnessi, Fig. 3) are in conflict with the pattern suggested by the combined nuclear DNA data. Therefore, we regard the evolutionary relationships among Oryctolagus, Caprolagus, Bunolagus, and Pentalagus as equivocal.

We anticipated that most of the mtDNA homoplasic characters were on the interior branches of the tree and the mtDNA data were contributing significantly to resolving the terminal associations among taxa. To test this we constructed topologies employing a Bayesian approach using a subset of taxa (including only the four Pronolagus, six Sylvilagus, twoNesolagus, and six Lepus species). These analyses resulted in 24 mtDNA trees within the 95% credibility interval, while 16 trees were found for the nuclear DNA data using the same taxa. Moreover, greater resolution resulted when all seven genes were combined (four unique trees were found).

Leporid Evolution

The data presented here suggest that the current leporid distribution is the result of at least nine dispersal events since the origin of the group in either North America or Asia. Three of these involved intercontinental exchanges between North America and Asia via the Bering strait, one between Asia and Europe, and at least three more between Europe/Asia and Africa. The evolution of the group is also characterized by at least two more recent dispersal events into South America (one by Lepus and another by Sylvilagus). Based on the molecular clock, it seems as if most of these events happened during the Miocene (between 14 and 8 MYA), the period during which all the modern leporid genera originated.

It has been suggested on palaeontological grounds that the Leporidae evolved either in Asia (Hibbard, 1963; Dawson, 1981) or North America (Voorhies and Timperley, 1997). Given the current distribution of the Leporidae, the DiVa analyses concur, suggesting that the ancestor of all modern leporids were either Asian or Asian and North American (C or CA; Fig. 4). It thus seems reasonable to suggest that the first intercontinental exchange giving rise to the present distribution of the leporids took place between North America and Asia. The molecular clock estimations suggest that this predates 14.00 (± 1.47) MYA. Although the direction of the first exchange is not clear from our analyses it is noteworthy that tectonic events during the middle Miocene epoch (∼ 20 MYA) resulted in the formation of land bridge connections between Eurasia and North America (Kummel, 1970). These geological changes were also accompanied by oscillations in Earth's orbit, resulting in severe global climate changes that included the melting of the Antarctic ice sheets (Zachos et al., 2001). Of particular importance to the evolution of the leporids is that the global ice volume was thought to have remained relatively stable and low from 26 to 15 MYA, thus isolating Asia from North America (Zachos et al., 2001). We hypothesize that the opening up of the Bering straight during the middle Miocene (approx 15 MYA) allowed for one leporid lineage to either established itself in North America (descendants of which include Romerolagus, Lepus, Sylvilagus, Brachylagus, Pentalagus, Caprolagus, Bunolagus, and Oryctolagus) or, alternatively in Asia (descendants of which include Nesolagus, Pronolagus, and Poelagus). Interestingly, it has been suggested that the ancestral leporid originated in forested environments (Voorhies and Timperley, 1997), and it is noteworthy that the most basal taxon within each clade (Nesolagus and Romerolagus) are confined to dense vegetation occurring at higher altitudes (Nowak, 1999).

Following this initial continental exchange between Asia and North America, the Asian ancestor extended its radiation into Africa giving rise to the extant Poelagus (in central Africa) and to Pronolagus (east and southern Africa; node E, Fig. 3). TheDiVa analyses indicated a vicariance event subsequent to this African immigration at approximately 11.3 MYA (± 1.53 MYA), separating the Asian Nesolagus from the African Poelagus/Pronolagus assemblage.

On the other hand, within North America the leporid ancestor gave rise to two lineages which diverged approximately 12.77 MYA (± 1.34 MYA): one giving rise to Romerolagus and the other to the lineage currently comprising the extant Lepus, Sylvilagus, Brachylagus, Pentalagus, Caprolagus, Bunolagus, and Oryctolagus. Shortly hereafter, between 11.8 MYA (± 1.24) MYA and 10.28 MYA (± 1.34 MYA), a second intercontinental exchange occurred between North America and Asia and this exchange was followed by a vicariance event leaving the ancestor of Sylvilagus and Brachylagus in North America and causing the establishment of ancestral Pentalagus, Caprolagus, Bunolagus, and Oryctolagus in Asia. The “cottontail lineage” (Sylvilagus and Brachylagus, node G) initially remained in North America (as initially did Lepus, see below) from where it colonized South America more recently. In this regard, it is noteworthy that our data place the pygmy rabbit, Brachylagus idahoensis, at the root of this clade, clearly distinct from the cottontails (Sylvilagus), giving credence to its separate generic status.

The Asian/European/African diversification of ancestral Pentalagus, Caprolagus, Bunolagus, and Oryctolagus happened more or less contemporaneously through a process of two dispersal and two vicariance events (Fig. 4). The uncertainties surrounding the exact branching pattern of this clade makes the biogeographic interpretation moot. Although the nuclear DNA topology made more geographic sense regarding the current distribution of the Japanese Pentalagus (it was also the only taxon that had a different placement when the two gene trees were compared—in the supermatrix analyses Pentalagus clustered with the geographically distant South African Bunolagus and in the combined nuclear phylogeny Pentalagus was placed basal to the geographically close Caprolagus from Asia), DiVa analyses based on the nuclear topology and the supermatrix respectively was equally parsimonious (reflecting two dispersals and two vicariance events). Irrespective which topology is correct, it is noteworthy that the isolation of Pentalagus from the remainder of this group is estimated at around 9.44 MYA (± 1.15 MYA) and this coincides closely with the Pliocene isolation of the Ryukyu Archipelago from the Asian continent (Yamada et al., 2002); the archipelago includes the Amami Oshima and Tokuno-shima islands that currently support Pentalagus. The reasons for the isolation of the remaining three genera, Caprolagus (India), Oryctolagus (Europe), and Bunolagus (South Africa) are unclear, but does suggest that up to now Africa was colonized twice, once before 11.3 MYA by the lineage that gave rise to the extant Poelagus and Pronolagus, and then more recently by the derived Bunolagus. Both Caprolagus and Bunolagus are now endangered, their distribution is restricted to a few remnant fragments of previously more widely distributed early successional riverine grassland habitats on the northern Indian subcontinent and in the South Africa karoo, respectively. Also, Pentalagus and Oryctolagus both share the unusual leporid trait of being colonial burrowers.

The third and final wave of dispersal out of North America into Asia was far more recent and was represented by a single genus, Lepus.Hibbard (1963) suggested that “the genus Lepus have arisen from both pro-Sylvilagus and Oryctolagus instead of from just one of the stocks.” Our molecular phylogeny is congruent with such a statement and strongly suggests that the divergence of the genus predates the development of the Sylvilagus/Oryctolagus clade (∼ 10.28 MYA; Fig. 3). Our molecular clock calibration places the divergence of the true hares, Lepus, at approximately 11.8 MYA, with four continental dispersal events probably occurring more recently. It has been suggested that the Lepus radiation coincided with the opening up the temperate grasslands, which was facilitated by the development of precocial young and cursorial adaptations suitable for an open habitat (enlargement of body size and strong hind legs; Corbet, 1986; Yamada et al., 2002). Although it has been argued that some Lepus species invaded North America more recently, as a result a secondary interchange among continents (Halanych et al., 1999), it is our thesis that the initial trigger for the continental dispersal of the genus Lepus was the global development of grasslands at ∼ 7 to 5 MYA and the formation of the west Antarctic ice sheet (6.5 MYA; Zachos et al., 2001). These factors once again promoted the development of land bridges that allowed dispersal from North America through Asia into Africa.

Acknowledgements

This work was supported by the University of Stellenbosch and the South African National Research Foundation (GUN 2053662). Wilbur Harrison is thanked for some laboratory assistance and Chris Ray, Fred Elder, Paula Jenkins (British Museum of Natural History), Ute Kryger, and Masashi Harada provided tissue for DNA extractions. Jack Sullivan, Luis Ruedas, and one anonymous reveiwer are thanked for valuable comments on an earlier draft of the manuscript and Jeffrey Thorne is thanked for assistance with the multidivergence software

References

Allard
M. W.
Honeycutt
R. L.
Nucleotide sequence variation in the mitochondrial 12S rRNA gene and the polyphyly of African mole-rats (Rodentia: Bathyergidae)
Mol. Biol. Evol.
 , 
1992
, vol. 
9
 (pg. 
27
-
40
)
Angerman
R.
Flux
J. E. C.
Chapman
J. A.
Smit
A. T.
Chapman
J. A.
Flux
J. E. C.
Lagomorph Classification
Rabbits, hares and pikas: Status conservation action plan
 , 
1990
Gland, Switzerland
International Union for Conservation of Nature and Natural Resources
(pg. 
7
-
13
Pages
Baker
R. H.
DeSalle
R.
Multiple sources of character information and the phylogeny of Hawaiian drosophilids
Syst. Biol.
 , 
1997
, vol. 
46
 (pg. 
654
-
673
)
Bollback
J. P.
Bayesian model adequacy and choice in phylogenetics
Mol. Biol. Evol.
 , 
2002
, vol. 
19
 (pg. 
1171
-
1180
)
Buckley
T. R.
Model misspecification and probabilistic tests of topology: Evidence from empirical data sets
Syst. Biol.
 , 
2002
, vol. 
51
 (pg. 
509
-
523
)
Buckley
T. R.
Arensburger
P.
Simon
C.
Chambers
G. K.
Combined data, Bayesian phylogenetics, and the origin of the New Zealand cicada genera
Syst. Biol.
 , 
2002
, vol. 
51
 (pg. 
4
-
18
)
Buckley
T. R.
Simon
C.
Chambers
G. K.
Exploring among-site rate variation models in a maximum likelihood framework using empirical data: Effects of model assumptions on estimates of topology, branch lengths, and bootstrap support
Syst. Biol.
 , 
2001
, vol. 
50
 (pg. 
67
-
86
)
Can
D. N.
Abramov
A. V.
Tikhonov
A. N.
Averianov
A. O.
Annamite striped rabbit Nesolagus timminsi in Vietnam
Acta Theriol.
 , 
2001
, vol. 
46
 (pg. 
437
-
440
)
Chapman
J. A.
Ceballos
G.
Chapman
J. A.
Flux
J. E. C.
The Cottontails. Pages 95–110 in Rabbits, hares and pikas: Status conservation action plan
 , 
1990
Gland, Switzerland
International Union for Conservation of Nature and Natural Resources
Chapman
J. A.
Cramer
K. L.
Dippenaar
N. J.
Robinson
T. J.
Systematics and biogeography of the New England cottontail, Sylvilagus transitionalis (Bangs, 1895), with the description of a new species from the Appalachian Mountains
Proc. Biol. Soc. Washington
 , 
1992
, vol. 
105
 (pg. 
841
-
866
)
Chapman
J. A.
Flux
J. E. C.
Chapman
J. A.
Flux
J. E. C.
Introduction and overview of the Lagomorpha. Pages 1–6 in Rabbits, hares and pikas: Status conservation action plan
1990
Gland, Switzerland
International Union for Conservation of Nature and Natural Resources
Cognato
A. I.
Vogler
A. P.
Exploring data interaction and nucleotide alignment in a multiple gene analysis of Ips (Coleoptera: Scolytinae)
Syst. Biol.
 , 
2001
, vol. 
50
 (pg. 
758
-
780
)
Corbet
G. B.
A review of classification in the family Leporidae
Acta Zool. Fennica
 , 
1983
, vol. 
174
 (pg. 
11
-
15
)
Dawson
M.
Later Tertiary Leporidae of North America. Univ
Kansas Paleont. Contr., Vertebrate
 , 
1958
, vol. 
6
 (pg. 
1
-
75
)
Dawson
M. R.
Myers
K.
MacInnes
C. D.
Evolution of modern leporids
Proceedings of the World Lagomorph Conference
 , 
1981
Ontario
University of Guelph Press
(pg. 
1
-
8
Pages
Diersing
V. E.
Anderson
S.
Knox Jones
J.
Jr.
Lagomorphs
Orders and Families of recent mammals of the world
 , 
1984
New York
John Wiley and Sons
(pg. 
241
-
254
Pages
Dowton
M.
Austin
A. D.
Increased congruence does not necessarily indicate increased phylogenetic accuracy-the behavior of the incongruence length difference test in mixed-model analyses
Syst. Biol.
 , 
2002
, vol. 
51
 (pg. 
19
-
31
)
Ellerman
J. R.
S. Morrison-Scott
T. C.
Checklist of palearctic and Indian mammals 1758–1946
 , 
1951
London
British Museum (Natural History)
Farris
J. S.
Kallersjo
M.
Kluge
A. G.
Bult
C.
Constructing a significance test for incongruence
Syst. Biol.
 , 
1995
, vol. 
44
 (pg. 
570
-
572
)
Flux
J. E. C.
Angermann
R.
Chapman
J. A.
Flux
J. E. C.
The hares and jackrabbits
Rabbits, hares and pikas: Status conservation action plan
 , 
1990
Gland, Switzerland
International Union for Conservation of Nature and Natural Resources
(pg. 
61
-
94
Pages
Frey
J. K.
Fisher
R. D.
Ruedas
L. A.
Identification and restriction of the type locality of the Manzano Mountains cottontail, Sylvilagus cognatus Nelson, 1907 (Mammalia: Lagomorpha: Leporidae)
Proc. Biol. Soc. Washington
 , 
1997
, vol. 
110
 (pg. 
329
-
331
)
Gatesy
J.
Amato
G.
Vrba
E.
Schaller
G.
DeSalle
R.
A cladistic analysis of mitochondrial ribosomal DNA from the Bovidae
Mol. Phylogenet. Evol.
 , 
1997
, vol. 
7
 (pg. 
303
-
319
)
Gatesy
J.
Matthee
C.
DeSalle
R.
Hayashi
C.
Resolution of a supertree/supermatrix paradox
Syst. Biol.
 , 
2002
, vol. 
51
 (pg. 
652
-
664
)
Gatesy
J.
Milinkovitch
M.
Waddell
V.
Stanhope
M.
Stability of cladistic relationships between Cetacea and higher-level artiodactyl taxa
Syst. Biol.
 , 
1999
, vol. 
48
 (pg. 
6
-
20
)
Gissi
C.
Gullberg
A.
Arnason
U.
The complete mitochondrial DNA sequence of the rabbit, Oryctolagus cuniculus
Genomics
 , 
1998
, vol. 
50
 (pg. 
161
-
169
)
Halanych
K. M.
Robinson
T. J.
Phylogenetic relationships of cottontails (Sylvilagus, Lagomorpha): Congruence of 12S rDNA and cytogenetic data
Mol. Phylogenet. Evol.
 , 
1997
, vol. 
7
 (pg. 
294
-
302
)
Halanych
K. M.
Robinson
T. J.
Multiple substitutions affect the phylogenetic utility of cytochrome b and 12S rDNA data: Examining a rapid radiation in leporid (Lagomorpha) evolution
J. Mol. Evol.
 , 
1999
, vol. 
48
 (pg. 
369
-
379
)
Halanych
K. M.
Demboski
J. R.
Jansen van Vuuren
B.
Klein
D. R.
Cook
J. A.
Cytochrome b phylogeny of North American hares and jackrabbits (Lepus, Lagomorpha) and the effects of saturation in outgroup taxa
Mol. Phylogenet. Evol.
 , 
1999
, vol. 
11
 (pg. 
213
-
221
)
Hibbard
C. W.
The origin of the P3 pattern of Sylvilagus, Caprolagus, Oryctolagus and Lepus
J. Mammal.
 , 
1963
, vol. 
44
 (pg. 
1
-
15
)
Hillis
D. M.
Mable
B. K.
Larson
A.
Davis
S. K.
Zimmer
E. A.
Hillis
D. M.
Moritz
C.
Mable
B. K.
Nucleic acids IV: Sequencing and cloning
Molecular systematics
 , 
1996
2nd ed.
Sunderland, MA
Sinauer Associates
(pg. 
321
-
381
Pages
Hillis
D. M.
Pollock
D. D.
McGuire
J. A.
Zwickl
D. J.
Is sparse taxon sampling a problem for phylogenetic inference?
Syst. Biol.
 , 
2003
, vol. 
52
 (pg. 
124
-
126
)
Huelsenbeck
J. P.
Testing a covariotide model of DNA substitution
Mol. Biol. Evol.
 , 
2002
, vol. 
19
 (pg. 
698
-
707
)
Huelsenbeck
J. P.
Ronquist
F.
MRBAYES: Bayesian inference of phylogenetic trees
Bioinformatics
 , 
2001
, vol. 
17
 (pg. 
754
-
755
)
Irwin
D. M.
Kocher
T. D.
Wilson
A. C.
Evolution of the cytochrome b gene of mammals
J. Mol. Evol.
 , 
1991
, vol. 
32
 (pg. 
128
-
144
)
Källersjö
M.
Albert
V. A.
Farris
J. S.
Homoplasy increases phylogenetic structure
Cladistics
 , 
1999
, vol. 
15
 (pg. 
91
-
93
)
Kelsey
C. R.
Crandall
K. A.
voevodin
A. F.
Different models, different trees: The geographic origin of PTLV-I
Mol. Phylogenet. Evol.
 , 
1999
, vol. 
13
 (pg. 
336
-
347
)
Kishino
H.
Thorne
J. L.
Bruno
W. J.
Performance of a divergence time estimation method under a probabilistic model of rate evolution
Mol. Biol. Evol.
 , 
2001
, vol. 
18
 (pg. 
352
-
361
)
Kocher
T. D.
Thomas
W. K.
Meyer
A.
Edwards
S. V.
Pääbo
S.
Villablanca
F. X.
Wilson
A. C.
Dynamics of mitochondrial DNA evolution in animals: Amplification and sequencing with conserved primers
Proc. Natl. Acad. Sci. U. S. A.
 , 
1989
, vol. 
86
 (pg. 
6196
-
6200
)
Kummel
B
History of the earth: An introduction to historical geology
 , 
1970
San Francisco
W. H. Freeman and Company
Matthee
C. A.
Burzlaff
J. D.
Taylor
J. F.
Davis
S. K.
Mining the mammalian genome for artiodactyl systematics
Syst. Biol.
 , 
2001
, vol. 
50
 (pg. 
367
-
390
)
Matthee
C. A.
Davis
S. K.
Molecular insights into the evolution of the family Bovidae: A nuclear DNA perspective
Mol. Biol. Evol.
 , 
2001
, vol. 
18
 (pg. 
1220
-
1230
)
Matthee
C. A.
Robinson
T. J.
Mitochondrial DNA differentiation among geographical populations of Pronolagus rupestris, Smith's red rock rabbit (Mammalia: Lagomorpha)
Heredity
 , 
1996
, vol. 
76
 (pg. 
514
-
523
)
Matthee
C. A.
Robinson
T. J.
Cytochrome b phylogeny of the family bovidae: Resolution within the alcelaphini, antilopini, neotragini, and tragelaphini
Mol. Phylogenet. Evol.
 , 
1999
, vol. 
12
 (pg. 
31
-
46
)
Meng
J.
Wyss
A. R.
Dawson
M. R.
Renjie
Z.
Primitive fossil rodent from Inner Mongolia and its implications for mammalian phylogeny
Nature
 , 
1994
, vol. 
370
 (pg. 
134
-
136
)
Morrone
J. J.
Crisci
J. V.
Historical biogeography: Introduction to methods
Ann. Rev. Ecol. Syst.
 , 
1995
, vol. 
26
 (pg. 
373
-
401
)
Murphy
W. J.
Eizirik
E.
Johnson
W. E.
Zhang
Y. P.
Ryder
O. A.
O'Brien
S. J.
Molecular phylogenetics and the origins of placental mammals
Nature
 , 
2001
, vol. 
409
 (pg. 
614
-
618
)
Murrell
A.
Campbell
N. J.
Barker
S. C.
A total-evidence phylogeny of ticks provides insights into the evolution of life cycles and biogeography
Mol. Phylogenet. Evol.
 , 
2001
, vol. 
21
 (pg. 
244
-
258
)
Nowak
R. M.
Walker's mammals of the worl
 , 
1999
6th ed.
Baltimore
John Hopkins University Press
Pääbo
S.
Wilson
A. C.
Polymerase chain reaction reveals cloning artifacts
Nature
 , 
1988
, vol. 
334
 (pg. 
387
-
388
)
Pollock
D. D.
Zwickl
D. J.
McGuire
J. A.
Hillis
D. M.
Increased taxon sampling is advantageous for phylogenetic inference
Syst. Biol.
 , 
2002
, vol. 
51
 (pg. 
664
-
671
)
Posada
D.
Crandall
K. A.
Modeltest: Testing the model of DNA substitution
Bioinformatics
 , 
1998
, vol. 
14
 (pg. 
817
-
818
)
Posada
D.
Crandall
K. A.
Selecting the best-fit model of nucleotide substitution
Syst. Biol.
 , 
2001
, vol. 
50
 (pg. 
580
-
601
)
Robinson-Rechavi
M.
Huchon
D.
RRTree: Relative-rate tests on a tree
Bioinformatics
 , 
2000
, vol. 
16
 (pg. 
296
-
297
)
Robinson
T. J.
Elder
F. F. B.
Chapman
J. A.
Evolution of chromosomal variation in cottontails, genus Sylvilagus (Mammalia: Lagomorpha). 2. Sylvilagus audubonii, S. idahoensis, S. nuttalli and S . palastris Cytogenet
Cell Genet.
 , 
1984
, vol. 
38
 (pg. 
282
-
289
)
Robinson
T. J.
Yang
F.
Harrison
W. R.
Chromosome painting refines the history of genome evolution in hares and rabbits (order Lagomorpha). Cytogenet
Genome Res.
 , 
2002
, vol. 
96
 (pg. 
223
-
227
)
Ronquist
F.
Dispersal-vicariance analysis: A new approach to the quantification of historical biogeography
Syst. Biol.
 , 
1997
, vol. 
46
 (pg. 
195
-
203
)
Ronquist
F.
Huelsenbeck
J. P.
MrBayes 3 Bayesian phylogenetic inference under mixed models
Bioinformatics
 , 
2003
, vol. 
19
 (pg. 
1572
-
1574
)
Springer
M. S.
Hollar
L. J.
Burk
A.
Compensatory substitutions and the evolution of the mitochondrial 12S rRNA gene in mammals
Mol. Biol. Evol.
 , 
1995
, vol. 
12
 (pg. 
1138
-
1150
)
Springer
M. S.
Murphy
W. J.
Eizirik
E.
O'Brien
S. J.
Placental mammal diversification and the Cretaceous-Tertiary boundary
Proc. Natl. Acad. Sci. U. S. A.
 , 
2003
, vol. 
100
 (pg. 
1056
-
1061
)
Su
C.
Nei
M.
Fifty-million year old polymorphism at an immunoglobulin variable region gene locus in the rabbit evolutionary lineage
Proc. Natl. Acad. Sci. U. S. A.
 , 
1999
, vol. 
96
 (pg. 
9710
-
9715
)
Sullivan
J.
Swofford
D. L.
Should we use model-based methods for phylogenetic inference when we know that assumptions about among-site rate variation and nucleotide substitution pattern are violated?
Syst. Biol.
 , 
2001
, vol. 
50
 (pg. 
723
-
729
)
Surridge
A. K.
Timmins
R. J.
Hewitt
G. M.
Bell
D. J.
Striped rabbits in Southeast Asia
Nature
 , 
1999
, vol. 
400
 pg. 
726
 
Swofford
D. L.
PAUP: Phylogenetic Analysis using Parsimony (and other methods). Version 4
 , 
2000
Sunderland, MA
Sinauer Associates
Thorne
J. L.
Kishino
H.
Divergence time and evolutionary rate estimation with multilocus data
Syst. Biol.
 , 
2002
, vol. 
51
 (pg. 
689
-
702
)
Thorne
J. L.
Kishino
H.
Painter
I. S.
Estimating the rate of evolution of the rate of molecular evolution
Mol. Biol. Evol.
 , 
1998
, vol. 
15
 (pg. 
1647
-
1657
)
Voorhies
M. R.
Timperley
C. L.
A new Pronotolagus (Lagomorpha: Leporidae) and other leporids from the Valentine railway quarriesw (Barstovian, Nebraska), and the archaeolagine-leporine transition
J. Ver. Paleon.
 , 
1997
, vol. 
17
 (pg. 
725
-
737
)
Waddell
P. J.
Cao
Y.
Hasegawa
M.
Mindell
D. P.
Assessing the Cretaceous superordinal divergence times within birds and placental mammals by using whole mitochondrial protein sequences and an extended statistical framework
Syst. Biol.
 , 
1999
, vol. 
48
 (pg. 
119
-
137
)
Whelan
S.
Lio
P.
Goldman
N.
Molecular phylogenetics: State-of-the-art methods for looking into the past
Trends Genet.
 , 
2001
, vol. 
7
 (pg. 
262
-
272
)
White
J. A.
North American Leporidae (Mammalia, Lagomorpha) from Late Miocene (Clarendonian) to late Pliocene (Blancan)
J. Ver. Paleon.
 , 
1991
, vol. 
1
 (pg. 
67
-
89
)
Whiteford
C. E. M.
Molecular phylogeny of the genus Pronolagus (Mammalia: Lagomorpha) and the use of morphological and molecular characters in the delineation of P. rupestris
 , 
1995
University of Pretoria
 
MSc Thesis
Yamada
F.
Takaki
M.
Suzuki
H.
Molecular phylogeny of Japanese Leporidae, the Amami rabbit Pentalagus furnessi, the Japanese hare Lepus brachyurus, and the mountain hare Lepus timidus, inferred from mitochondrial DNA sequences
Genes Genet. Syst.
 , 
2002
, vol. 
77
 (pg. 
107
-
116
)
Yang
Z.
PAML: A program package for phylogenetic analysis by maximum likelihood
Comput. Appl. Biosci.
 , 
1997
, vol. 
13
 (pg. 
555
-
556
)
Yang
Z.
Yoder
A. D.
Comparison of Likelihood and Bayesian methods for estimating divergence times using multiple gene loci and calibration points, with application to a radiation of cute-looking mouse lemus species
Syst. Biol.
 , 
2003
, vol. 
52
 (pg. 
1
-
12
)
Yu
N.
Zheng
C.
Zhang
Y. P.
Li
W. H.
Molecular systematics of pikas (genus Ochotona) inferred from mitochondrial DNA sequences
Mol. Phylogenet. Evol.
 , 
2000
, vol. 
16
 (pg. 
85
-
95
)
Zachos
J.
Pagani
M.
Sloan
L.
Thomas
E.
Billups
K.
Trends, rhythms, and aberrations in global climate 65 Ma to present
Science
 , 
2001
, vol. 
292
 (pg. 
686
-
693
)
Zwickl
D. J.
Hillis
D. M.
Increased taxon sampling greatly reduces phylogenetic error
Syst. Biol.
 , 
2002
, vol. 
51
 (pg. 
588
-
598
)