The Volcano Rabbit in the Phylogenetic Network of Lagomorphs.

The order Lagomorpha unifies pikas (Ochotonidae) and the hares plus rabbits (Leporidae). Phylogenetic reconstructions of the species within Leporidae based on traditional morphological or molecular sequence data provide support for conflicting hypotheses. The retroposon presence/absence patterns analyzed in this study revealed strong support for the broadly accepted splitting of lagomorphs into ochotonids and leporids with Pronolagus as the first divergence in the leporid tree. Furthermore, the retroposon presence/absence patterns nested the rare volcano rabbit, Romerolagus diazi, within an unresolved network of deeper leporid relationships and provide the first homoplasy-free image of incomplete lineage sorting and/or ancestral hybridization/introgression in rapidly radiated Leporidae. At the same time, the strongest retroposon presence/absence signal supports the volcano rabbit as a separate branch between the Pronolagus junction and a unified cluster of the remaining leporids.

The volcano rabbits (Romerolagus diazi), also known as zacatuche in their Mexican homeland, are one of the smallest rabbits with no visible tail and short round ears (Fa and Bell 1990). Contrary to cosmopolitan hares (the genus Lepus), they populate only the central part of the Mexican transvolcanic belt at elevations between 2,800 and 4,250 m (Cervantes et al. 1990; Fa and Bell 1990). This habitat probably acted as a refuge for volcano rabbits, allowing them to survive during the Pleistocene (Ceballos et al. 2010). Today, however, this strong habitat specialization makes them vulnerable to environmental perturbations. With just a little more than 2,000 remaining, their status is classified as endangered in the IUCN Red List of Threatened Species (https://www.iucnredlist.org/details/19742/0; last accessed December 17, 2018). Destruction of their narrow natural habitats and illegal hunting heightens this species’ vulnerability to extinction.

The volcano rabbit represents a monotypic genus within the family Leporidae (rabbits and hares), which includes 11 extant genera and, together with Ochotonidae (pikas), forms the order Lagomorpha. The origin of modern Leporidae lies in North America, where they differentiated in the early and middle Miocene. From North America, leporids dispersed to northern Asia across Beringia ∼8 Ma, or possibly earlier, and spread worldwide during the late Miocene (Flynn et al. 2014). However, there is still a lack of consensus concerning the phylogenetic relationships among leporid genera and the phylogenetic position of the volcano rabbit in the leporid tree (fig. 1). Although it was suggested that the volcano rabbit has a number of ancestral morphological traits (Velazquez et al. 1993), recent morphological analyses of lagomorph skulls did not provide significant phylogenetic resolution due to convergent evolution (Ge et al. 2015; Kraatz and Sherratt 2016). Based on molecular data, only the basal position of the rock rabbit (Pronolagus) is assured (mitochondrial data, Halanych and Robinson 1999; Ge et al. 2013; nuclear and nuclear+mt data, Matthee et al. 2004; supertree, Rolland et al. 2014). In some studies, Pronolagus is grouped with Nesolagus (Matthee et al. 2004; Ge et al. 2013) and Poelagus (Matthee et al. 2004), or alternatively with Pentalagus (Rolland et al. 2014). The volcano rabbit was proposed to form the second basal branch of leporids (supermatrix data including five nuclear and two mitochondrial gene fragments; Matthee et al. 2004). This analysis revealed the divergence time of the volcano rabbit from other leporids to be ∼12.8 Ma, which is close to estimates based on 12S rDNA (Halanych and Robinson 1999). The subsequent supermatrix data of Robinson and Matthee (2005) assembled Matthee’s sequence data and insertion/deletion patterns, together with morphological, cytogenetic, and geographical features resulted in the same tree topology as that provided by Matthee et al. (2004). A supertree built by Bininda-Emonds et al. (2007), updated by Fritz et al. (2009), and further refined using the Kuhn et al. polytomy resolver (Kuhn et al. 2011; Rolland et al. 2014) also provided support for the volcano rabbit diverging from the leporid root following the basal split of the Pronolagus branch. The karyotype of the volcano rabbit is similar to hares (chromosome number of 2n = 48 and G-banding pattern) but differs from most other rabbits and was proposed to be an ancestral feature that was retained in the ancient volcano rabbit lineage (van der Loo et al. 1979; Robinson et al. 1981). However, this consensus notwithstanding, a recent study based on three mitochondrial genes revealed a considerably more terminal position for the volcano rabbit, placing it as a sister group to a Bunolagus/Pentalagus assemblage, after Oryctolagus and Sylvilagus separated (Ge et al. 2013). Moreover, Ge et al. estimated the volcano rabbit divergence to be more recent (∼5.7 Ma) than posited in previous studies.

. 1.

—Contradictory conclusions of the phylogenetic position of Romerolagus in the lagomorph tree. Left: Phylogenetic tree of Lagomorpha derived from nuclear and mitochondrial sequences (Matthee et al. 2004); Right: An alternative proposal of the phylogenetic position of Romerolagus derived from mitochondrial sequences (Ge et al. 2013). Only lagomorph genera that were investigated in the present study are presented on the trees.

Irrespective of these considerations, however, it is believed that the leporid lineages diversified over a short time (Halanych and Robinson 1999), and consequently the discordant phylogenetic position of the volcano rabbit may result not only from homoplasy affecting data sets but also by hemiplasy occurring as a result of incomplete lineage sorting (ILS). Hemiplasy is a phenomenon that is reflected in gene tree–species tree discordance. It arises from the phylogenetic sorting of ancestral polymorphisms over successive speciation nodes, followed by subsequent random fixation, or loss, that results in homoplasy-like consequences for lineages, despite the fact that character states are genuinely homologous (see definition of hemiplasy in Avise and Robinson [2008]; Robinson et al. [2008]). An alternative but not exclusive process possibly underlying the conflicting phylogenetic signals is ancestral hybridization, which may result in introgression or hybrid speciation (e.g., Baack and Rieseberg 2007).

Most phylogenetic marker systems cannot easily distinguish ILS/ancestral hybridization from other random noise (homoplasy); however, the retroposon marker system presents a more reliable exception (Doronina et al. 2019). Shared retroposon insertions at orthologous loci in the genomes of two species indicate a diagnostic integration in their common ancestral genome and thus their close phylogenetic relationship. As with any phylogenetic marker system, retroposon presence/absence signals that were polymorphic in ancestral populations at the time of speciation could, potentially, lead to the occurrence of hemiplasious signals contradicting the species-tree; furthermore, their distribution among lineages can be affected by ancestral hybridization. However, carefully analyzed retroposon presence/absence patterns are well suited to recognizing such ILS/hybridization-produced signals because of their virtually homoplasy-free nature and may therefore be especially useful in evaluating phylogenetic relationships that are blurred by ILS and/or hybridization effects (Kuritzin et al. 2016). Moreover, and particularly pertinent to the present study, lagomorph-specific C families of short interspersed elements (CSINEs) were present and active across the full diversification of this order (Kriegs et al. 2010), further underscoring the possibility that they may be useful in resolving the unsettled phylogenetic position of the volcano rabbit within Leporidae. The Kriegs et al. (2010) study was hitherto the only one in which a retroposon marker system was used to investigate lagomorph phylogenetic relationships; however, the volcano rabbit was not included in their sampling.

As there are currently only two lagomorph genomes that have been sequenced (the American pika, Ochotona princeps, and the European rabbit, Oryctolagus cuniculus), we performed two in silico screenings for retroposon insertions: 1) present in both Oryctolagus and Ochotona and 2) present in Oryctolagus but absent in Ochotona. We took the advantage of the available assembled and annotated genome of Oryctolagus, while in the previous lagomorph retroposon study only low-quality Oryctolagus trace sequences were available (Kriegs et al. 2010). This enabled us to perform an exhaustive, genome-wide screening of the Oryctolagus genome and to use Oryctolagus annotation data to increase the number and length of extracted introns (up to 1,100 vs. <800 nt based on the human assembly in Kriegs et al. [2010]). We extracted 17,831 introns of Oryctolagus, 536 of which contained potential phylogenetically informative CSINE insertions (compared with just 74 intron loci available in Kriegs et al. [2010]). Fifty-eight loci contained CSINE1 elements, 384 loci CSINE2s, and 94 loci CSINE3s. As CSINE2s represent the oldest of the CSINEs and were generally active at the root of Lagomorpha and Leporidae (Kriegs et al. 2010), they were less informative for younger splits, and so we have now focused mainly on CSINE1 and CSINE3 elements that were active during Leporidae diversification. After constructing manual alignments of these loci and conducting experimental identification of retroposon presence or absence states in lagomorph species lacking genome sequence data, we collected those that were amplifiable in a representative set of species including the volcano rabbit. We identified 23 novel diagnostic retroposon markers from 19 loci (fig. 2 and table 1) that were not previously analyzed by Kriegs et al. (2010). To better identify the phylogenetic position of the volcano rabbit, we also reanalyzed three markers (current markers 6, 7, and 19) that were earlier found by Kriegs et al. (2010; see their supplementary table S1, markers 10, 9, 8, respectively) on the branch leading to Oryctolagus, Bunolagus, Sylvilagus, and Lepus (see their fig. 2), but whose presence/absence patterns were not previously determined for the volcano rabbit (fig. 2 and table 1).

. 2.

—(A) Retroposon-based phylogenetic tree of Lagomorpha. Gray balls represent markers found in this study and those found by Kriegs et al. (2010) and reanalyzed in the present study regarding their Romerolagus presence/absence states (table 1); open circles are markers taken from Kriegs et al. (2010) without reanalysis. Markers representing incongruent presence/absence patterns are shown to the right of the figure. For two of the four markers placing Romerolagus as the second basal group, presence/absence states in Bunolagus were not identified. P-values from the one-directional KKSC insertion significance test are provided under the nodes that were significantly supported by retroposon presence/absence data. (B) Neighbor-net analysis (SplitsTree) of generic-level phylogenetic relationships within Leporidae based on retroposon presence/absence patterns. Lepus sp. combines all clear presence/absence states of the five investigated Lepus species. Bootstrap values (numbers) provided only for groups supported by retroposon markers (gray balls).

Table 1

Presence/Absence Patterns of 26 Informative Retroposon Markers Detected in This Study

Marker	Element	Lagomorpha											Out
		Opr	Pcr	Rdi	Ocu	Bmo	Sfl	Lca	Leu	Lto	Lsa	Lcap
1b	CSINE2	1	1	1	1	1	1	?	1	1	1	?	0*
2	CSINE1	0	1	1	1	1	1	1	1	1	?	1	0
3	CSINE2	0	1	1	1	?	1	1	1	1	1	1	0
4b	CSINE1	0	1	1	1	1	1	1	1	1	1	1	0
5	CSINE2	0	1	1	1	1	1	1	?	?	?	?	d
6(K)	CSINE3	0	0	1	1	1	1	1	1	1	1	1	0
7(K)	CSINE1	0	0	1	1	1	1	1	1	1	1	1	0
8	CSINE1	0	0	1	1	1	1	1	1	1	1	1	0
9b	CSINE1	0	0	1	1	1	1	1	1	1	1	1	0
10a	CSINE3	0	0	1	1	1	1	1	1	1	1	1	0
11	CSINE3	0	0	1	1	1	1	1	1	1	1	1	0
12	CSINE1	0	0	1	1	1	1	1	1	1	1	1	0
13	CSINE3	0	0	1	1	?	1	1	1	?	?	?	0*
1a	CSINE1	0	0	1	1	1	1	?	1	1	1	?	0*
10b	CSINE3	0	0	0	0	0	0	1	1	1	1	1	0
4a	CSINE3	0	0	0	0	0	0	1	1	1	1	1	0
14	CSINE3	0	0	0	0	0	0	d	1	1	1	1	0
9a	CSINE3	0	0	0	0	0	0	0	0	0	1	1	0
15	CSINE3	0	0	1	1	0	0	0	0	0	?	?	0
16	CSINE1	0	?	1	1	1	0	0	?	?	0	0	0
17	CSINE3	0	0	0	1	1	1	0	0	0	0	0	0
18	CSINE3	0	0	0	1	1	1	0	0	0	0	0	0
19(K)	CSINE3	0	0	0	1	1	1	1	1	?	1	1	0
20	CSINE3	0	0	0	1	?	1	?	1	1	1	1	0
21	CSINE3	0	0	0	1	1	1	1	1	?	?	1	0*
22	CSINE3	0	0	0	1	?	1	1	?	1	1	1	0

Note.—(1) indicates the presence of retroposed element, (0) indicates its absence, (?) denotes lack of sequence information, (d) denotes deletions in the retroposon insertion region. Markers denoted by letters are located in one locus (e.g., 1a and 1b in locus 1). Opr, Ochotona princeps; Pcr, Pronolagus crassicaudatus; Rdi, Romerolagus diazi; Ocu, Oryctolagus cuniculus; Bmo, Bunolagus monticularis; Sfl, Sylvilagus floridanus; Lca, Lepus californicus; Leu, Lepus europaeus; Lto, Lepus townsendii; Lsa, Lepus saxatilis; Lcap, Lepus capensis; Out, outgroup. In the column Out (0) indicates absence in Cavia porcellus, (0*)—in Mus musculus. The markers 6, 7, and 19 marked as (K) were initially found by Kriegs et al. (2010) (their markers 10, 9, and 8, respectively) and reanalyzed in the present study regarding the Romerolagus presence/absence state.

In addition to the markers identified by Kriegs et al. (2010) (fig. 2, open circles), we found one new CSINE2 marker that supports the monophyly of the Lagomorpha (i.e., pikas and leporids) and four new markers (two CSINE1s and two CSINE2s) supporting leporid monophyly (fig. 2, gray balls; P < 0.004, [8 0] KKSC insertion significance test; Kuritzin et al. 2016). These data also confirm and refine the Kriegs et al. (2010) finding that CSINE2 elements comprise the oldest lagomorph SINEs by showing that their retroposition activity appears to have ceased prior to the first leporid diversification (table 1).

Within Leporidae, the basal position of Pronolagus as a sister group to all other leporids is significantly supported by nine retroposon markers (fig. 2P < 0.003, [9 0] KKSC insertion significance test). Four diagnostic retroposon markers support the position of Romerolagus as the second diverged leporid group; however, we were unable to verify the presence/absence state of two of them in Bunolagus (fig. 2). We note that the only available leporid genome of Oryctolagus restricted us to perform Oryctolagus-directed screenings and did not allow us to analyze alternative phylogenetic affiliations of Romerolagus (e.g., Romerolagus–Bunolagus sister group relationship). However, we found clear indications for a more basal position of Romerolagus rather than a terminal position (fig. 2). Although the χ2 test was not significant for the basal position of Romerolagus (P > 0.5), a neighbor-net analysis yielded a bootstrap support of 98.6%, indicating strong support for Romerolagus as the second basal branch (fig. 2B), which is in agreement with the tree topology of Matthee et al. (2004) rather than that of Ge et al. (2013). Moreover, we also found two markers supporting the Oryctolagus+Bunolagus+Sylvilagus group that does not contradict the second basal position of Romerolagus and provides an additional argument in favor of Matthee et al. (2004) over Ge et al. (2013) in the phylogenetic tree reconstruction.

However, we also found one conflicting signal merging Romerolagus with Oryctolagus and one marker supporting the Romerolagus+Oryctolagus+Bunolagus group (fig. 2). Interestingly, the previous retroposon phylogeny of lagomorphs did not detect any contradictory markers (Kriegs et al. 2010), which might be explained by the smaller number of loci analyzed and a restricted species sampling excluding the volcano rabbit. Given the rapid radiation in leporids (Halanych and Robinson 1999) and taking into account the virtually homoplasy-free nature of retroposons (Doronina et al. 2019), the incongruent phylogenetic signals found in our study are most probably due to ILS, which was shown to accompany the diversifications of many mammalian groups (Churakov et al. 2009; Nishihara et al. 2009; Doronina et al. 2015; Feigin et al. 2018). Alternatively, the conflicting signals might also be due to ancestral hybridization that may have led to introgression or hybrid speciation and may have played an essential role in lagomorph species evolution. Previously, interspecific mitochondrial and nuclear introgression events were detected in Lepus (Liu et al. 2011; Melo-Ferreira et al. 2014). The presence of ancestral ILS or ancestral hybridization in leporid diversifications is a very important finding for future investigations of leporid relationships and the phylogenetic position of Romerolagus. Given that complete genomes comprise a mosaic of different blocks that may have independent histories and provide contradicting phylogenetic signals (Pääbo 2003; Doronina et al. 2017), further resolution of the leporid phylogeny will depend on the addition of genome assemblies for the remaining representatives of this important mammalian group.

Materials and Methods

Based on the European rabbit (Or.cuniculus, oryCun2) genome annotation, we extracted sequences of all short introns (300–1,100 nt) and 50 nt for each of their conserved flanking exons for primer design using the UCSC Table browser (http://genome.ucsc.edu/cgi-bin/hgTables; last accessed December 17, 2018). Following a local RepeatMasker run (http://www.repeatmasker.org; last accessed December 17, 2018), we extracted introns containing full-length CSINE elements. We used the National Center for Biotechnology Information (NCBI) BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi; last accessed December 17, 2018) and UCSC Genome Browser Database BLAT (http://genome.ucsc.edu/cgi-bin/hgBlat; last accessed December 17, 2018) to search these sequences for orthologous loci in the genome of the American pika (O.princeps) and outgroups (guinea pig [Cavia porcellus] or mouse [Mus musculus] in case guinea pig was not available). We constructed manual alignments and chose potentially informative loci containing the following presence(+)/absence(−) retroposon patterns: rabbit(+)/pika(+)/outgroup(−) and rabbit(+)/pika(−)/outgroup(−). We identified presence/absence patterns of the insertions in leporids experimentally using the same species that were included by Kriegs et al. (2010) along with Romerolagus diazi and Lepus townsendii (table 1). We designed conserved polymerase chain reaction primers (supplementary file S1, Supplementary Material online) and performed polymerase chain reaction amplification, cloning, and sequencing using standard protocols (Farwick et al. 2006). In addition to markers found in the present screens, we also included in our analyses three markers from Kriegs et al. (2010) that were potentially informative for the Romerolagus phylogenetic position and experimentally identified their presence/absence states in Romerolagus (table 1, loci 6, 7, and 19). All alignments are provided as supplementary file S2, Supplementary Material online. This project has been deposited at GenBank under the accession MK078357–MK078509.

To examine the significance of the data supporting the phylogenetic nodes, we used the KKSC insertion significance test (Kuritzin et al. 2016) located at http://retrogenomics.uni-muenster.de:3838/KKSC_significance_test/; last accessed December 17, 2018. Because the sequence of only one leporid genome (Or. cuniculus) is currently available, we were restricted to performing screens and statistical analyses of markers only in one direction. In addition, we performed a χ2 test according to Waddell et al. (2001). To reconstruct the leporid phylogenetic network, we also performed neighbor-net analysis of a presence/absence (1/0) data matrix (supplementary file S3, Supplementary Material online) in SplitsTree (Huson and Bryant 2006; version 4.13.1) using the uncorrected P-distance and default settings (ignoring ambiguous states; using normalize option; allowing a 100% missing data level per site). In the matrix, [1] denotes a presence state of a retroposon, [0] denotes the absence, and [?] denotes lack of sequence information or a deletion in the retroposon insertion region. Only markers representing leporid intergeneric relationships were included in this analysis.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.