Illuminating the obscured phylogenetic radiation of South American Sylvilagus Gray, 1867 (Lagomorpha: Leporidae)

Abstract

Largely shallow and putatively explosive divergences in the family Leporidae (rabbits and hares; order: Lagomorpha) have resulted in phylogenetic relationships that remain currently unresolved. These rapid radiations in different branches of the leporid tree have resulted in conflicting phylogenetic hypotheses. However, this phylogenetic incongruence may also result from inadequate taxon or character sampling, due to the high number of extinct and difficult to sample extant species, and highly conserved morphological characters. Sylvilagus (cottontail rabbits) constitute about 30% of the known extant leporid species. New species are routinely being recognized, and phylogenetic relationships with respect to other leporid genera, and within the genus, have failed to be recovered with certainty. Within Sylvilagus, the South American S. brasiliensis is the most widespread and poorly known taxon, likely comprising multiple species. Here, we reanalyze previously published molecular data from phylogenetic studies on Leporidae, focusing on the S. brasiliensis group, and assessing phylogenetic relationships using bifurcating trees and split networks to identify phylogenetic regions with polytomies. We estimate differentiation and phylogenetic relationships of molecular lineages within the S. brasiliensis group. Our analyses suggest that this group contains a number of divergent taxa, well differentiated from other cottontail species. We discern two major polytomies during leporid diversification. The first, at the base of the leporid radiation, likely resulted from a combination of hard (rapid radiation) and soft polytomies (high number of unsampled extinct species). The second polytomy likely resulted from a rapid radiation during the initial diversification of the genus Sylvilagus. We conclude that only a molecular phylogeny based on a broader taxonomic representation will fully resolve leporid phylogeny.

Divergências recentes, e provavelmente explosivas, na família de mamíferos Leporidae (coelhos e lebres; ordem: Lagomorpha) resultaram em relações filogenéticas ainda por resolver. Rápidas radiações em diferentes ramos da árvore dos leporídeos têm originado hipóteses filogenéticas conflitantes. Contudo, a incongruência filogenética pode ser resultado de amostragem taxonômica ou de caracteres inadequadas, devido a um alto número de taxa extintos e de espécies viventes difíceis de amostrar, ou a caracteres morfológicos bastante conservados. O gênero Sylvilagus (coelhos americanos) compreende cerca de 30% das espécies de leporídeos viventes. Novas espécies de Sylvilagus estão continuamente sendo descritas e, nem as relações filogenéticas com outros gêneros de leporídeos, nem as relações dentro do próprio gênero foram ainda recuperadas com certeza. Dentro de Sylvilagus, o táxon sul americano S. brasiliensis é o mais amplamente distribuído e menos conhecido, provavelmente incluindo múltiplas espécies. Neste trabalho, examinamos hipóteses filogenéticas moleculares previamente publicadas para Leporidae, expandindo a representação taxonômica do grupo S. brasiliensis e verificando as relações filogenéticas usando árvores bifurcadas e redes filogenéticas para identificar regiões de politomia. Nós estimamos a diferenciação e relações filogenéticas das linhagens moleculares dentro do grupo S. brasiliensis. As nossas análises sugerem que o grupo de espécies S. brasiliensis contém taxa bem divergentes das outras espécies do gênero. Além disso, observamos duas politomias principais ocorridas durante a diversificação dos leporídeos. A primeira, na base da sua radiação, provavelmente resultou da combinação de rápida radiação e elevado número de espécies extintas não amostradas. A segunda politomia provavelmente provém da rápida radiação durante a fase inicial de diversificação de Sylvilagus. Nós concluímos que apenas filogenias moleculares baseadas numa representação taxonômica mais ampla conseguirão resolver a filogenia dos leporídeos.

molecular phylogenetics, Neotropical diversity, rapid radiations, supernetworks, Tapeti

Assessment of relationships and phylogenetic reconstructions of pikas, hares, and rabbits (Lagomorpha) have historically been difficult. These mammals were described as—and considered—rodents for the first century and a half of the Linnean era (Gidley 1912; Arnason et al. 2002; Asher 2005; Kraatz et al. 2009; Ruedas et al. 2018). Similarly, the shallow phylogenetic relationships among hares and rabbits (Leporidae) also have been inconsistently resolved (Yamada et al. 2002; Matthee et al. 2004; Robinson and Matthee 2005; Ruf 2014). Not even combining cytogenetic, molecular, and morphological data has shed significant light on their relationships, because homoplasy, convergence, and autapomorphies have, in concert, hampered attaining more conclusive phylogenetic inferences (Robinson and Matthee 2005; Ruf 2014; Ge et al. 2015).

Introgression (Melo-Ferreira et al. 2009; Liu et al. 2011), natural selection and adaptation (Carneiro et al. 2012; Melo-Ferreira et al. 2014), and rapid radiations in different regions of the leporid tree as currently resolved (Halanych and Robinson 1999; Halanych et al. 1999; Robinson and Matthee 2005; Ruedas et al. 2017) have resulted in complex evolutionary histories and distinct phylogenetic hypotheses. The phylogenetic incongruencies could also have resulted from soft polytomies (i.e., insufficient data—Coddington and Scharff 1996). The most comprehensive molecular phylogenetic matrices of Lagomorpha published to date either used a multilocus approach and 25 of at least 67 extant leporid species (Matthee et al. 2004: figure 3: molecular supermatrix; Burgin et al. 2018), or were based on mtDNA only, comprising 31 species (Ge et al. 2013: figure 3). Resolving soft polytomies in the phylogeny of Leporidae is a difficult task, not only because of the number of extant species and the combination of a conservative bauplan with potentially homoplasious morphological characters, but also due to the elusive behavioral habits and low densities of the majority of the species, which hinders their sampling and study, as well as because of the large number of extinct species. Up to 30 genera of extinct leporids have been described to date (Flynn et al. 2014). Thus, both patterns of phylogenetic relationships and species diversity within the family, extinct and extant, remain unresolved (Jin et al. 2010; Pelletier et al. 2015; Ruedas 2017).

Sylvilagus (Gray, 1867), the American cottontail rabbits, occurs from southeastern Canada to northern Argentina. This genus is of particular interest with respect to the aforementioned topics: phylogenetic relationships in particular, both among congeners as well as other leporid genera, have been inconsistently determined (Halanych and Robinson 1997; Matthee et al. 2004; Ge et al. 2013; Bonvicino et al. 2015; Ruedas et al. 2017). A rapid speciation rate within the genus was observed by Ruedas et al. (2017) to have been accompanied with low extinction rates. Recognized Sylvilagus species constitute 22% of all Lagomorpha and 33% of the extant leporids’ diversity (Burgin et al. 2018; Smith et al. 2018). However, taxonomic uncertainties remain, with new species having been recently recognized (Ruedas 1998, 2017; Durant and Guevara 2001), and the number of IUCN Data Deficient and Threatened species in Sylvilagus being the highest among lagomorphs (Verde Arregoitia et al. 2015).

Until recently, Sylvilagus brasiliensis (Linnaeus, 1758), the South American cottontail, was considered the most widespread American lagomorph, occurring from Mexico to Argentina (Silva Jr. et al. 2005; AMCELA et al. 2008). However, ongoing detailed analyses are discriminating among taxa in what is now recognized as a complex of species (Ruedas and Salazar-Bravo 2007; Ruedas 2017; Ruedas et al. 2017). Previously, over 25 subspecies have, at one time or another, been subsumed within S. brasiliensis, making this one of the most polytypic taxa in the genus (Hoffmann and Smith 2005). Morphological and genetic data, allied with biogeographic and ecological niche-based analyses, suggest that S. brasiliensis (sensu stricto) is restricted to a narrow forested area within the Pernambuco Endemism Center of the northeastern Atlantic Forest biome, Brazil (Ruedas et al. 2017). Sylvilagus tapetillus (Thomas, 1913), formerly included within S. brasiliensis, is either extinct or restricted to a minute coastal plain region (Serra do Mar ecoregion—da Silva et al. 2004; Dinerstein et al. 2017), in the southeastern Atlantic Forest biome (Bonvicino et al. 2015; Ruedas et al. 2017). Specimens from central and southern Atlantic Forest (formerly assigned either to S. brasiliensis or S. tapetillus), and Cerrado taxa (putatively S. brasiliensis, S. minensis [Thomas, 1901], or S. paraguensis [Thomas, 1901]) currently have an unclear taxonomic status (Ruedas et al. 2017). Occurrence records of putative S. brasiliensis specimens within the Amazon forest and Caatinga were not considered in the IUCN Red List, nor have these been discriminated with a subspecific epithet (Silva Jr. et al. 2005; AMCELA et al. 2008; Dantas et al. 2016). In addition, populations putatively ascribed to S. brasiliensis are considered Endangered in Rio Grande do Sul state, southern Brazil (state decree 51797, 8 September 2014). Moreover, the European brown hare (Lepus europaeus) and European rabbit (Oryctolagus cuniculus) both are expanding their ranges across South America (Novillo and Ojeda 2008; Bonino et al. 2010; Costa and Fernandes 2012; Galende 2014). Thus, cottontail rabbits within the region are most likely already competing for resources with these larger-sized competitors, given that both are highly successful invasive species and constitute major conservation challenges wherever they have been introduced (Lees and Bell 2008; Novillo and Ojeda 2008).

Traditional representation of phylogenetic hypotheses as bifurcating trees might be uninformative or even misleading when evolutionary histories include rapid radiations (Suh 2016), such as those hypothesized for leporids (Halanych and Robinson 1999; Halanych et al. 1999; Robinson and Matthee 2005; Ruedas et al. 2017). Instead, comparing bifurcating trees to phylogenetic networks can help identify problematic nodes and possible explosive radiation events, thereby ruling out soft polytomies (Huson and Bryant 2006; Suh et al. 2015; Hahn and Nakhleh 2016). Supernetworks are particularly useful as they enable combining data sets with different terminal taxa, resulting in a more comprehensive set of phylogenetic evidence (Huson et al. 2004).

Here, we reanalyze previously published molecular data from phylogenetic studies on leporids (Matthee et al. 2004; Ge et al. 2013; Ruedas et al. 2017), focusing on the S. brasiliensis species group, and assess phylogenetic relationships based on bifurcating trees and split networks. Our purpose was to identify otherwise cryptic divergent molecular lineages previously ascribed to S. brasiliensis, test the hypothesis of a common origin of these lineages, and describe evolutionary relationships among lineages (putative species), their congeners, and other genera of leporids.

Materials and Methods

Taxon sampling.

We sampled the recently described (Ruedas et al. 2017) neotype of S. brasiliensis (Universidade Federal de Pernambuco [UFPE] no. 1740), 11 Brazilian specimens from the S. brasiliensis group never before included in any previous phylogenetic inference, and a L. europaeus specimen (Table 1; Fig. 1). We also used sequences retrieved from the GenBank from 81 leporids used in other studies (Appendix I) to perform the phylogenetic analyses detailed below. Newly collected S. brasiliensis specimens were obtained from different phytophysiognomic regions and latitudes, and could potentially be assigned to species distinct from S. brasiliensis (sensu Ruedas et al. 2017; Table 1). These specimens were collected under the auspices of Brazilian permit number 43504 (MMA/ICMBio/SISBIO) and following ASM guidelines (Sikes et al. 2016).

Fig. 1.

The most accepted distribution of Sylvilagus brasiliensis (sensu lato; shaded area—AMCELA et al. 2008). Point locations for all samples used from the S. brasiliensis group (except for S. gabbi and S. dicei), and Lepus europaeus (L). New samples correspond to localities 1–12, and samples from previous studies correspond to point locations 13–22 (see Table 1 and Appendix I for more information). Symbols correspond to the clades in Table 2 (black triangle, SbrNeot: S. brasiliensis Neotype UFPE 1740; black and gray diamonds, Smn, S. minensis; black circles, EAM, Eastern Amazon; black squares, SW, Southwest; black cross, CR, Costa Rica; and the four divergent specimens, respectively, named in the map).

Table 1.

Geographic sampling information. Taxonomic assignments were attributed following Ruedas et al. (2017) and this study. A sample from Lepus europaeus (*) was also included. The ecoregion assigned was the main type in the area of occurrence, following Dinerstein et al. (2017). Numbers correspond to point locations depicted in Fig. 1. Institutional abbreviations: UFPE, Federal University of Pernambuco; UFMG, Federal University of Minas Gerais; UERJ, Rio de Janeiro State University; LABEC, Laboratory of Evolutionary Biology and Conservation of Vertebrates, São Paulo University; DEPAVE, Department of Parks and Green Areas of São Paulo—Faunal Division; LBM, Laboratory of Molecular Biology; MPEG, Museu Paraense Emílio Goeldi; MCNU, Museum of Natural Sciences—Lutheran University of Brazil. Neot.: Neotype. Ecoregion abbreviations: PIE, Pernambuco interior forests; Ce, Cerrado; BIF, Bahia interior forests; SMCF, Serra do Mar coastal forests; APAF, Alto Paraná Atlantic forest; MGTDF, Mato Grosso tropical dry forests, TPMF, Tocantins/Pindare moist forest; **collection locality imprecise: this region is in a transition zone between TPMF, Ce, and Maranhão Babaçu forests.

	Sample	Institution	Species		Locality	Ecoregion
			Ruedas et al. (2017)	This study
1	UFPE 1740	UFPE	S. brasiliensis Neot.		Mata da Privativa, Paudalho, Pernambuco	PIE
2	M1380	UFMG	S. minensis	S. minensis 1	Parque das Mangabeiras, Belo Horizonte, Minas Gerais	Ce
3	M1778	UFMG	S. minensis	S. minensis 1	São Gonçalo do Rio Abaixo, Minas Gerais	BIF
4	M1770	UFMG	S. minensis	S. minensis 1	Nova Xavantina, Mato Grosso	Ce
5	LG479	UERJ	S. minensis	S. minensis 1	Jequitinhonha Valley, Bahia	BIF
6	DPV 53580	LABEC/DEPAVE	S. minensis	S. minensis 2	São Paulo state	SMCF
7	SP01	LBM/MPEG	S. minensis	S. minensis 2	Parque Anhanguera, São Paulo, São Paulo	SMCF
8	RS01	LBM/MPEG	S. minensis	S. paraguensis	Parque Estadual do Turvo, Rio Grande do Sul	APAF
9	M1796	UFMG	Sylvilagus sp.	Sylvilagus sp. 1	Paraupebas, Pará	MGTDF
10	TOC013	LABEC/DEPAVE	Sylvilagus sp.	Sylvilagus sp. 2	Maranhão state	**
11	MPEG 45455	MPEG	Sylvilagus sp.	Sylvilagus sp. 2	Benjamin Constant, Pará	TPMF
12	MPEG 45456	MPEG	Sylvilagus sp.	Sylvilagus sp. 2	Mirasselvas, Pará	TPMF
L	MCNU 1561*	MCNU	Lepus europaeus		São José das Missões, Rio Grande do Sul

Table 2.

Within and between clades uncorrected p-distances for mitochondrial 12S gene in individuals of the Sylvilagus brasiliensis group. Standard error estimates are shown above the diagonal. SbrNeot: S. brasiliensis Neotype UFPE 1740; Smn, S. minensis; EAM, Eastern Amazon; SW, Southwest; CR, Costa Rica. A * symbol indicates that the group is represented by a single sample. These groups are represented by different symbols and colors in the map of Fig. 1.

	SbrNeot*	Smn	EAM	SW	M1796*	MN24041	ROM105515*	CR	CR1hsr*
SbrNeot*		0.002	0.002	0.003	0.003	0.006	0.004	0.006	0.007
Smn	0.003	0.000	0.004	0.003	0.003	0.006	0.003	0.006	0.007
EAM	0.004	0.008	0.002	0.004	0.004	0.007	0.006	0.007	0.008
SW	0.010	0.010	0.010	0.004	0.003	0.006	0.004	0.006	0.007
M1796*	0.007	0.007	0.010	0.010		0.006	0.004	0.005	0.007
MN24041*	0.025	0.023	0.027	0.027	0.024		0.005	0.003	0.006
ROM105515*	0.011	0.008	0.015	0.014	0.010	0.023		0.005	0.007
CR	0.026	0.024	0.027	0.028	0.024	0.006	0.023	0.003	0.006
CR1hsr*	0.042	0.040	0.040	0.043	0.041	0.027	0.038	0.027

Laboratory procedures.

We extracted genomic DNA from fresh tissue samples and samples from tissues preserved in ethanol using the DNeasy Blood & Tissue kit (Qiagen, Valencia, California). We amplified the molecular markers used by Matthee et al. (2004): two mitochondrial genes (cytochrome b [Cytb] and 12S ribosomal RNA [12S]) and five nuclear introns (mast cell growth factor [MGF], protein kinase C iota [PRKCI], thyroid-stimulating hormone beta [THY], beta spectrin nonerythrocytic protein 1 [SPTBN], and thyroglobin [TG]). PCR conditions and primers followed Matthee et al. (2004) (with minor modifications; Supplementary Data SD1) and Ruedas et al. (2017) for mtDNA loci. PCR products were visually inspected after electrophoresis in a 1% agarose gel. Positive results were purified using PEG8000 2.5M (Hawkins et al. 1994). Sequences for both strands were obtained in an ABI 3130 automatic sequencer, using the Big Dye Terminator v3.01 kit, following the manufacturer’s protocol (Applied Biosystems, Foster City, California).

Genetic diversity patterns within the S. brasiliensis group.

All new sequences generated for this study were deposited in GenBank, with accession numbers MH115201–MH115262. New sequences were first visually inspected in Chromas (http://technelysium.com.au/wp/chromas/) and aligned using ClustalW (Thompson et al. 1994) as implemented in Bioedit 7 (Hall 2011). The PHASE algorithm implemented in DnaSP 5 was used to describe haplotypes for nuclear molecular markers (Stephens and Scheet 2005; Librado and Rozas 2009). MCMC chains were run for 1,000 iterations with a burn-in of 100. Genetic diversity indices and Tajima’s D neutrality tests were obtained using DnaSP (Librado and Rozas 2009). For the complete S. brasiliensis group, we used the 12S locus only (Appendix I) to estimate mean uncorrected p-distances in MEGA 5, running 1,000 bootstrap replications (Tamura et al. 2011).

Data sets.

Different sets of species and loci were used to obtain improved resolution of problematic nodes and to assess the effect of soft polytomies on leporid phylogeny, as follows: Data set 1—expanded molecular supermatrix: included all loci, and only Matthee et al. (2004) sequences and our new Sylvilagus sequences (n = 38 ingroup taxa; 4,792 base pairs [bp]); Data set 2—Cytb data set (n = 56 ingroup taxa; 582 bp); Data set 3—12S data set (n = 72 ingroup taxa; 742 bp); and Data set 4—Cytb and 12S data sets combined (n = 93 ingroup taxa; 1,324 bp; Supplementary Data SD1). Cytochrome b is the locus with greatest taxonomic representation in Leporidae, with 47 species represented in GenBank, whereas 12S is the most comprehensive data set for Sylvilagus species, with 39 species.

Phylogenetic analyses.

We used sequences of Tamiasciurus hudsonicus (Rodentia, Sciuridae) and Ochotona princeps (Lagomorpha, Ochotonidae) as outgroups in our phylogenetic analyses (Appendix I; Matthee et al. 2004). Sequences were aligned as above. When necessary, sequences were trimmed to decrease missing data. Maximum likelihood (ML) trees were obtained using RAxML 8.1.21 (Stamatakis 2014) with the graphical interface raxmlGUI 1.3.1 (Silvestro and Michalak 2012), generating 1,000 replicates with the rapid bootstrap algorithm (options -f a), and branch lengths calculated independently for each partition. Best substitution models implemented in RAxML were compared, and the best set of partitions was estimated in PartitionFinder (Lanfear et al. 2012) for each data set. The Cytb data set was divided into two partitions (codon positions 1 + 2 and 3), and Cytb alignment in the expanded molecular supermatrix data set was divided by codon position. All other alignments corresponded to a single partition. A GTR-GAMMA mutation model was used in both ML and Bayesian inferences (BI). Multilocus coalescent Bayesian trees were reconstructed using BEAST 1.8 (Drummond et al. 2012) at the CIPRES Science Gateway 3.3 (www.phylo.org—Miller et al. 2010). Initial runs were carried out under a relaxed-molecular-clock approach with an uncorrelated lognormal distribution and a speciation Yule process as tree prior (Yule 1925; Drummond et al. 2012; Heled and Drummond 2015). Final runs consisted of at least two independent chains with a length of 10⁸ generations. Effective sample sizes (ESS) and the convergence of the chains were assessed using Tracer 1.5; the initial 10% of each independent run was excluded as burn-in (Drummond and Rambaut 2007). Trees were combined with LogCombiner (Drummond et al. 2012), and tree topologies assessed using TreeAnnotator (Drummond et al. 2012) and FigTree 1.4 (http://tree.bio.ed.ac.uk/software/figtree/). File formats were converted using ALTER (http://sing.ei.uvigo.es/ALTER/—Glez-Peña et al. 2010).

Rapid diversification events in leporid evolution.

In order to visualize conflicting nodes, we followed the recommendations of Hahn and Nakhleh (2016) and Suh (2016), and constructed split networks for each molecular marker, and supernetworks for combined data sets (expanded molecular supermatrix, all sequences from all loci and taxa combined, separate Cytb and 12S sequences from taxa represented in the expanded molecular supermatrix, and full Cytb and 12S data sets combined; Appendix I). We used SplitsTree 4.14 (Huson et al. 2004; Huson and Bryant 2006), with 500 runs to estimate supernetworks.

Results

Genetic diversity patterns within the S. brasiliensis group.

We sequenced a total of 2,152 bp of mtDNA and 2,815 bp of nDNA introns. Several heterozygote sequences were found in all nDNA (except for SPTBN), and only two indels were observed, in 12S and PRKCI, respectively. No departures from neutrality were found. For the newly acquired S. brasiliensis species group sequences, the number of segregating sites varied from 4 (in PRKCI, SPTBN, and THY) to 87 (in Cytb), the number of haplotypes from 4 (SPTBN) to 12 (Cytb), and nucleotide diversity from 1.11 × 10⁻³ (THY) to 2.79 × 10⁻² (Cytb).

Pairwise differences between samples assigned to the S. brasiliensis group reached 5% (uncorrected p-distance; $\overline{x}$ = 0.016, range 0.00–0.048; Supplementary Data SD2). Sylvilagus dicei (Harris, 1932) and S. gabbi (Allen, 1877), previously assigned to S. brasiliensis (Diersing 1981; Ruedas and Salazar-Bravo 2007), were not included in p-distances calculations (but see Ruedas et al. 2017, appendix 2). The neotype of S. brasiliensis differed on average by 1.3% from the other samples of the group (0.003–0.042; Supplementary Data SD2). Lower distances were observed between South American clades and the neotype (⁠$\overline{x}$ < 1%), and higher distances were estimated between the neotype and the Central American clades (⁠$\overline{x}$ > 4%; Table 2).

In most phylogenies, neotype sequences clustered as sister to samples putatively assigned to S. minensis, except for RS01 (Fig. 2; Supplementary Data SD3, node A). This latter sample, from southern Atlantic Forest, Brazil, grouped with samples previously attributed to S. b. paraguensis and S. b. peruanus (Hershkovitz, 1950) for the 12S data set (Fig. 2C; Supplementary Data SD3, node E; SW group in Fig. 1 and Table 2). DPV 53580 and SP01, from southeastern Brazil, were recovered as monophyletic in two data sets (Figs. 2A and 2B; Supplementary Data SD3, node B). Specimens from the Amazonia formed two clades (Figs. 2B and 2C; Supplementary Data SD3, nodes C and D, corresponding to eastern Amazon EAM group and M1796 in Table 2, respectively). A Central American clade grouped samples from Costa Rica and MN24041 (Sbr24041) from Pernambuco (Fig. 2C; Supplementary Data SD3, node F), the latter identified as S. brasiliensis, and S. gabbi and S. dicei, as previously observed by Ruedas et al. (2017).

Fig. 2.

Phylogenetic trees for Leporidae species inferred from maximum likelihood for the extended molecular supermatrix (A), Cytb (B), and 12S (C) data sets. Values above branches indicate bootstrap support (≥ 50) for each node, and letters reference nodes and clades discussed in the text. Labels include specimen voucher if more than one sample was attributed to the same species; new samples have their locality of origin in brackets as in Table 1 and Fig. 1 (see also Appendix I for a full list of sequences used). Bayesian phylogenies are given in Supplementary Data SD3.

Evolutionary relationships among taxa in the S. brasiliensis group, their congeners, and the other leporids.

Bayesian inferences and ML trees presented considerable differences both in tree topologies and nodal support (Fig. 2; Supplementary Data SD3). BI for the expanded molecular supermatrix data set had little nodal support (Supplementary Data SD3). Relationships among genera were inconsistent across analyses. However, polytypic genera were supported as monophyletic with the exception of those represented in the expanded molecular supermatrix data set. Brachylagus idahoensis was supported as the sister clade to either Sylvilagus or Lepus, depending on the partitioning of the data and taxonomic representation. As mentioned above, samples assigned to the S. brasiliensis group were paraphyletic within the framework of the 12S data set (Fig. 2C; Supplementary Data SD3C). Central American samples grouped with S. dicei and S. gabbi (Fig. 2C; Supplementary Data SD3C, node F), while our novel samples grouped in an exclusively South American clade. The sister clade to the S. brasiliensis group likewise was inconsistently resolved. Either all other Sylvilagus, or S. floridanus (Allen, 1890) + S. obscurus (Chapman et al. 1992), alternatively had support as the sister clade under distinct phylogenetic analytical frameworks (Fig. 2; Supplementary Data SD3).

Rapid diversification events in leporid evolution.

Supernetworks for the expanded molecular supermatrix data set presented several zones of polytomy. We hypothesize that this may in part be due to missing data for the nDNA loci (Figs. 3A and 3B; Appendix I), given that using only mtDNA sequences from our data set decreased the number of polytomies (Fig. 3C). A supernetwork for Cytb and 12S data sets combined, comprising about twice the number of species represented in the expanded molecular supermatrix data set, resulted in only two polytomy zones (Fig. 3D). A polytomy zone is observed at the base of Leporidae, the other at the basal diversification of Sylvilagus.

Fig. 3.

Supernetworks for Leporidae based on (A) expanded molecular supermatrix, (B) all sequences combined, (C) Cytb and 12S from the expanded molecular supermatrix, and D) Cytb and 12S full data sets combined. In order to increase resolution, not all Sylvilagus samples were labeled. The origin of the new samples named is given in brackets as in Table 1 and Fig. 1. ▲ Neotype for S. brasiliensis (fully labeled networks are available upon request from the corresponding author).

Discussion

Genetic distances between what are considered valid species within Sylvilagus have been demonstrated to be variable. Recently diverged species (sister taxa) might be indistinguishable using mtDNA (e.g., S. obscurus versus S. transitionalis [Bangs, 1895]—Litvaitis et al. 1997; Ruedas et al. 2017), but others may differ by up to 5% (e.g., S. aquaticus [Bachman, 1837] versus S. palustris [Bachman, 1837]—Ruedas et al. 2017). As expected, we observed geographic partitioning: there was high divergence between Central and South American clades nominally constituted by putative S. brasiliensis (p-distances > 2.5%). Our analyses further revealed the existence of at least three divergent molecular lineages in Central America (Fig. 2, node F), more closely related to S. gabbi and S. dicei, both acknowledged as valid species-level taxa (Diersing 1981; Ruedas and Salazar-Bravo 2007). Our analyses also confirm the existence of multiple molecular lineages in the South American clade, among which the neotype stands out as a distinct lineage (Fig. 2; Supplementary Data SD3, node A), supporting the hypothesis that more comprehensive sampling remains to be undertaken. Nevertheless, samples attributed to S. minensis (sensu Ruedas et al. 2017; see Table 1) are closer to the neotype (p-distance = 0.03; Table 2), but appear to segregate from it, and into three distinct molecular lineages: one including different specimens from distinct ecoregions in central Brazil, another from the Serra do Mar coastal forests in the southeastern Atlantic Forest (DPV 53580 and SP01; Fig. 1, gray diamonds; Fig. 2; Supplementary Data SD2, node B), and a third from Alto Paraná Atlantic Forest in the southern Atlantic Forest (RS01). The first two clades should both be assigned to S. minensis, since they are not reciprocally monophyletic (Fig. 2; Supplementary Data SD3, clade A). RS01, which falls within a southwestern clade (Fig. 2; Supplementary Data SD3, node E), likely represents S. paraguensis (see also Ruedas et al. 2017: figures A20 and A24). Notwithstanding, several lineages remain to be excluded from the S. brasiliensis species group in this western South America clade (Ruedas 2017; Ruedas et al. 2017). Individuals from east and southeast Amazonia appear to constitute two distinct molecular lineages (Fig. 2; Supplementary Data SD3, nodes C and D). Table 1 summarizes our taxonomic conclusions on the new specimens analyzed.

Our analyses, particularly those using split networks (Fig. 3), demonstrate that the S. brasiliensis species assemblage, here including S. dicei and S. gabbi, represents a divergent evolutionary lineage, well differentiated from all other species in the genus. We also are confident in refuting the subgenus Tapeti (Gray, 1867) as a taxonomic entity, as we found no molecular evidence for grouping S. brasiliensis with the S. palustris + S. aquaticus lineage (see also Bonvicino et al. 2015; Ruedas et al. 2017). Furthermore, contrary to the most recent morphometric analyses focused on Central America and the north Andean region (Diersing and Wilson 2017), our results support the hypothesis of Ruedas et al. (2017) that S. andinus (Thomas, 1897) is a distinct species-level taxon that does not cluster with Central American S. brasiliensis specimens, nor within the S. brasiliensis group (Figs. 2C, 3B, and 3D; Supplementary Data SD3C).

Our supernetworks enable a more comprehensive phylogenetic hypothesis for Sylvilagus. These results show that the basal relationships among most cottontail species correspond to a zone of polytomy, even for data sets with a better taxonomic representation (Fig. 3), supporting the hypothesis of a rapid radiation in cottontails (Ruedas et al. 2017). The burst of diversification in Sylvilagus appears to have affected the current basal phylogenetic resolution of the genus, and this radiation has not slowed down (Ruedas et al. 2017). Moreover, our data corroborate previous findings that most cottontail species consistently cluster in pairs of sister species, such as S. palustris + S. aquaticus, S. floridanus + S. obscurus, and S. nuttallii (Bachman, 1837) + S. audubonii (Baird, 1858) (Matthee et al. 2004; Robinson and Matthee 2005; Ruedas et al. 2017). However, recovering phylogenetic relationships among these species pairs, as well as the sister species—if present—of the S. brasiliensis group has still not been achieved (this study; Bonvicino et al. 2015; Ruedas et al. 2017). Our data corroborate that the S. brasiliensis group is represented by far more than a pair of species (S. gabbi, S. dicei, S. brasiliensis, S. paraguensis, and others yet to be clarified—Ruedas et al. 2017). A better taxonomic understanding of polytypic species, including S. floridanus, S. bachmani (Waterhouse, 1839), and S. nuttallii, is required for a better understanding of the evolutionary history of the genus, as shown by the phylogenetic trees we obtained (Fig. 2; Supplementary Data SD3). Although morphological characters have overall low phylogenetic signal (Robinson and Matthee 2005; Ruf 2014; Ge et al. 2015), analyses of the third lower premolar and meticulous analysis of cranial characters have been successful in establishing species limits within Sylvilagus (Ruedas and Salazar-Bravo 2007) and appear to work well even in other closely related species (Ruedas 1998, 2017).

Brachylagus has been identified as the sister genus to Sylvilagus, albeit inconsistently (this study; Robinson and Matthee 2005; Ge et al. 2013; Ruedas et al. 2017). In this case, not even our approach using phylogenetic networks resolved this riddle. Yet the results of our supernetwork analyses show that Brachylagus and all the other leporid genera are included in a basal leporid polytomy (Fig. 3), justifying the support for monophyly in the different genera, despite the obscured phylogenetic relationships among them (Yamada et al. 2002; Matthee et al. 2004; Robinson and Matthee 2005; Ruf 2014). Therefore, on the one hand, the molecular variability assessed to date might be insufficient to resolve relationships among the Leporidae, particularly under the scenario of extremely rapid diversification that our results support (Halanych and Robinson 1999; Robinson and Matthee 2005; Ruedas et al. 2017): a rapid radiation alone could account for unresolved phylogenetic relationships (Suh 2016). On the other hand, increasing taxonomic representation by adding species in Sylvilagus and Lepus—i.e., decreasing soft polytomies—has had little effect in resolving the core zone of polytomy: the basal radiation of Leporidae. This might be due to the large number of leporids, from numerous genera, that became extinct during the Pleistocene (Ge et al. 2013; Flynn et al. 2014). Notwithstanding, the reduction of the polytomy basal to Lepus with more hare species sampled (Figs. 3C versus D) suggests that an increase in taxonomic representation of leporid taxa is favorable to decreasing the polytomies in their phylogeny, in spite of the scenario of rapid radiation.

In fact, for Lepus, supernetworks for 12S and Cytb data sets combined, the most comprehensive taxa data set for the genus, underscore the dissimilarity among recognized species (Fig. 3D), in particular, the dissimilarities between co-occurring or geographically proximal Lepus species. For instance, the relationships between L. americanus and the other North American and Mexican hares (Ramírez-Silva et al. 2010), or L. brachyurus from Japan and L. mandshuricus from the Asian mainland, support the existence of numerous bursts of radiation (Matthee et al. 2004; Ge et al. 2013). Moreover, and although bifurcating trees add limited or inconsistent support for most Lepus relationships, we unambiguously discriminated several clades in the networks, including monotypic lineages.

Our data have shed some light on the explosive basal radiation of Leporidae. The results of our analyses demonstrate the necessity to expand efforts at increased taxonomic and molecular sampling of additional leporid species. It is becoming increasingly evident that only by improving the taxonomic representation of species across the numerous branches of the leporid tree, it will be possible to clarify their phylogeny.

Nomenclature statement.

A life science identifier (LSID) number was obtained for this publication: urn:lsid:zoobank.org:pub:B727FE72-D2F8-4FE3-80D3-B2B5A2C7386D.

Supplementary Data

Supplementary data are available at Journal of Mammalogy online.

Supplementary Data SD1.—PCR conditions.

Supplementary Data SD2.—Pairwise uncorrected p-distances (lower diagonal) and standard error (upper diagonal) for mitochondrial 12S gene between samples attributed to the Sylvilagus brasiliensis group.

Supplementary Data SD3.—Phylogenetic trees for Leporidae species inferred from Bayesian inference for the extended molecular supermatrix (A), Cytb (B), and 12S (C) data sets. Values above branches indicate posterior probabilities (≥ 50) for each node, and letters reference nodes discussed in the text. Labels include specimen voucher if more than one sample was attributed to the same species, and new samples have their locality of origin in brackets as in Table 1 and Fig. 1 (see Appendix I for a full list of sequences used).

Acknowledgments

We thank D. A. de Moraes (Universidade Federal de Pernambuco), L. Geise (Universidade Estadual do Rio de Janeiro), J. Summa (Departamento de Parques e Áreas Verdes de São Paulo – Divisão de Fauna, DEPAVE), F. Santos and F. Perini (Universidade Federal de Minas Gerais), A. Cristoff (Universidade Luterana do Brasil do Rio Grande do Sul), and A. Bessa for access to specimens under their care; DEPAVE-SP staff (particularly M. Nardi and L. Lopes), L. Oliveira, U. Schulz, and M. Ferraz (Universidade do Vale do Rio dos Sinos), D. Meller, staff from Parque Estadual do Turvo (particularly V. Grutzmann), N. Simbera (Projeto Quatis), P. Avelino (UFMG), G. Dantas (Pontifícia Universidade Católica de Minas Gerais), A. Pavan (Universidade de São Paulo), and M. de Vivo (Museu de Zoologia da Universidade de São Paulo) for field support; and staff from MPEG for laboratory support. Fundação Amazônia de Amparo a Estudos e Pesquisas (ICAAF 009/2015) and Conselho Nacional de Desenvolvimento Científico e Tecnológico (MCTI/CNPq/Universal 447460/2014-5) funded field and lab work. SMS was funded by a fellowship from the Programa Nacional de Pós Doutorado/Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (PNPD/CAPES) at the Department of Zoology in the Museu Paraense Emílio Goeldi/Universidade Federal do Pará (PPGZOOL/MPEG/UFPA), and LHS by a scholarship from the Programa Institucional de Bolsas de Iniciação Científica (Conselho Nacional de Desenvolvimento Científico e Tecnológico) at the Department of Zoology MPEG/UFPA. LAR received support for portions of this work from National Science Foundation (grant DEB-0616305). We also acknowledge the comments and suggestions from the anonymous reviewers, which greatly improved this work.

Appendix I

Sample vouchers, GenBank accession numbers, and references for all sequences used in this study.

Voucher/isolate	Speciesa	CYTB	12S	MGF	PRKCI	TG	THY	SPTBN	References	Locality Fig. 1
UFPE1740	Sbr	MH115201	KU057257	MH115213	MH115222	MH115239	MH115251	MH115230	Ruedas et al. (2017); this study	1
M1380	Smn	MH115202	MH115190	MH115214		MH115240	MH115252		This study	2
M1770	Smn	MH115203	MH115191	MH115215	MH115223	MH115241	MH115253	MH115231	This study	4
M1778	Smn	MH115204	MH115192	MH115216	MH115224	MH115242	MH115254	MH115232	This study	3
LG479	Smn	MH115205	MH115193	MH115217	MH115225	MH115243	MH115255	MH115233	This study	5
DPV53580	Smn	MH115206	MH115194			MH115244		MH115234	This study	6
SP01	Smn	MH115207	MH115195			MH115245		MH115235	This study	7
RS01	Smn	MH115208	MH115196	MH115218	MH115226	MH115246	MH115256	MH115236	This study	8
M1796	Ssp	MH115209	MH115197	MH115219	MH115227	MH115247	MH115257	MH115237	This study	9
TOC013	Ssp	MH115210	MH115198			MH115248			This study	10
MPEG45455	Ssp	MH115211	MH115199	MH115220	MH115228	MH115249		MH115238	This study	11
MPEG45456	Ssp	MH115212	MH115200	MH115221	MH115229	MH115250	MH115258		This study	12
MCNU1561	Leu	MH115259	MH115260		MH115261	MH115262			This study	L
	Bid	AY292721	AY292695	AY292745	AY292772	AY292840	AY292669	AY292814	Matthee et al. (2004)
	Bmn	AY292718	AY292692	AY292742	AY292769	AY292837	AY292666	AY292811	Matthee et al. (2004)
	Chs	AY292719	AY292693	AY292743	AY292770	AY292838	AY292667	AY292812	Matthee et al. (2004)
	Lam	AY292733	AY292707	AY292757	AY292784	AY292852	AY292681	AY292826	Matthee et al. (2004)
	Lcf	AY292731	AY292705	AY292755	AY292782	AY292850	AY292679	AY292824	Matthee et al. (2004)
	Lcp	AY292732	AY292706	AY292756	AY292783	AY292851	AY292680	AY292825	Matthee et al. (2004)
	Lsx	AY292730	AY292704	AY292754	AY292781	AY292849	AY292678	AY292823	Matthee et al. (2004)
	Ltm	AY292728	AY292702	AY292752	AY292779	AY292847	AY292676	AY292821	Matthee et al. (2004)
	Ltw	AY292729	AY292703	AY292753	AY292780	AY292848	AY292677	AY292822	Matthee et al. (2004)
	Nnt		AY292709	AY292759	AY292786	AY292854	AY292683		Matthee et al. (2004)
	Ntm		AY292710	AY292760	AY292787	AY292855	AY292684	AY292828	Matthee et al. (2004)
	Ocn	AY292717	AY292691	AY292741	AY292768	AY292836	AY292665	AY292810	Matthee et al. (2004)
	Pcr	AY292738	AY292714	AY292764	AY292791	AY292859	AY292688	AY292832	Matthee et al. (2004)
	Pfr	AY292720	AY292694	AY292744	AY292771	AY292839	AY292668	AY292813	Matthee et al. (2004)
	Pmr				AY292767	AY292835			Matthee et al. (2004)
	Prn	AY292737	AY292713	AY292763	AY292790	AY292858	AY292687	AY292831	Matthee et al. (2004)
	Prp	AY292735	AY292711	AY292761	AY292788	AY292856	AY292685	AY292829	Matthee et al. (2004)
	Psa	AY292736	AY292712	AY292762	AY292789	AY292857	AY292686	AY292830	Matthee et al. (2004)
	Rdi	AY292734	AY292708	AY292758	AY292785	AY292853	AY292682	AY292827	Matthee et al. (2004)
	Saq	AY292726	AY292700	AY292750	AY292777	AY292845	AY292674	AY292819	Matthee et al. (2004)
	Sau	AY292722	AY292696	AY292746	AY292773	AY292841	AY292670	AY292815	Matthee et al. (2004)
	Sfl	AY292724	AY292698	AY292748	AY292775	AY292843	AY292672	AY292817	Matthee et al. (2004)
	Snt	AY292723	AY292697	AY292747	AY292774	AY292842	AY292671	AY292816	Matthee et al. (2004)
	Sob	AY292725	AY292699	AY292749	AY292776	AY292844	AY292673	AY292818	Matthee et al. (2004)
	Spl	AY292727	AY292701	AY292751	AY292778	AY292846	AY292675	AY292820	Matthee et al. (2004)
GK227	Lce	AB687531		LC131921		LC132018	LC132115	LC132212	Kinoshita et al. (2012)
s1	Lsn	HM233081		HM233721		HM233484			Liu et al. (2011)
ME	Lmn	HM233010		HM233344		HM233393			Liu and Zhang, Liu et al.
C1	Lcm	HM233066		HM233619		HM233377			Liu et al. (2011)
Lar3	Lar	JN037365			JN037012			JN037040	Melo-Ferreira et al. (2011)
Lcr1	Lci	JN037356			JN037003			JN037031	Melo-Ferreira et al. (2011)
Lot1	Lot	JN037366			JN037013			JN037041	Melo-Ferreira et al. (2011)
LbAKT243	Lbr	AB514237							Nunome et al. (2010)
DT	Lha	KF723323							Kong et al. (2016)
1	Lpg	KF723324							Kong et al. (2016)
O1	Loi	HM233064							Liu et al. (2011)
ASK 3518	Srb	HQ143448							Nalls et al. (2012)
	Lcs	AY942569							Melo-Ferreira et al. (2005)
NK 4504	Lal	AF010157							Halanych et al. (1999)
NK 3800	Lcl	AF010158							Halanych et al. (1999)
P151	Lgr	KJ943805							Lopes et al. (2014)
ECOSUR:2048	Lfl	KT308125							Álvarez-Castañeda and Lorenzo (2017)
RF-8	Lin	HM222713							Naidu et al. (2012)
	Ltb	LC073697							Shan and Liu (2016)
	Str	AF034256							Snyder and Husband
	Ltl	AY649617							Xiang, Yang and Xia
QCAZ10893	San		KU057258						Ruedas et al. (2017)
AK11178	Sau		KU057237						Ruedas et al. (2017)
MVZ206386	Sbc		KU057239						Ruedas et al. (2017)
MSB40683	Sbc		KU057238						Ruedas et al. (2017)
ROM105515	Sbr		KU057236						Ruedas et al. (2017)	13
MVZ153492	Sbr		KU057243						Ruedas et al. (2017)	14
MN24041	Sbr		KU057254						Ruedas et al. (2017)	15
MSB55948	Sbr		KU057240						Ruedas et al. (2017)	16
TK61307	Sbr		KU057242						Ruedas et al. (2017)	17
TTU79706	Sbr		KU057241						Ruedas et al. (2017)	18
EPN954419	Sbr		KU057228						Ruedas et al. (2017)	19
CR1hsr	Sbr		KU057219						Ruedas et al. (2017)	20
CRIV1	Sbr		KU057222						Ruedas et al. (2017)	21
CRIV2	Sbr		KU057223						Ruedas et al. (2017)	21
CRIV4	Sbr		KU057224						Ruedas et al. (2017)	21
CRIV5	Sbr		KU057225						Ruedas et al. (2017)	21
CRIV6	Sbr		KU057226						Ruedas et al. (2017)	21
EM1556	Sbr		KU057227						Ruedas et al. (2017)	22
TTU114374	Sdc		KU057256						Ruedas et al. (2017)
TK147518	Sdc		KU057251						Ruedas et al. (2017)
AK11511	Sfl		KU057246						Ruedas et al. (2017)
MVZ154373	Sfl		KU057231						Ruedas et al. (2017)
ASNHC_2330	Sfl		KU057218						Ruedas et al. (2017)
CR26gpv	Sfl		KU057220						Ruedas et al. (2017)
MSB89673	Sfl		KU057244						Ruedas et al. (2017)
hidra008	Sfl		KU057229						Ruedas et al. (2017)
IIBT349	Sfl		KU057253						Ruedas et al. (2017)
NP310	Sfl		KU057235						Ruedas et al. (2017)
MSB158806	Sgb		KU057232						Ruedas et al. (2017)
MSB158807	Sgb		KU057233						Ruedas et al. (2017)
PSU4944	Snt		KU057255						Ruedas et al. (2017)
AK11516	Sob		KU057248						Ruedas et al. (2017)
AK11529	Sob		KU057247						Ruedas et al. (2017)
PSU4960	Spl		KU057249						Ruedas et al. (2017)
AK11525	Str		KU057250						Ruedas et al. (2017)
	Opr	AY292716	AY292690	AY292740	AY292766	AY292834	AY292664	AY292809	Matthee et al. (2004)
	Tst	AY292715	AY292689	AY292739	AY292765	AY292833	AY292663	AY292808	Matthee et al. (2004)

^aBid, Brachylagus idahoensis; Bmn, Bunolagus monticularis; Chs, Caprolagus hispidus; Lal, Lepus alleni; Lam, L. americanus; Lar, L. arcticus; Lbr, L. brachyurus; Leu, L. europaues; Lcf, L. californicus; Lcl, L. callotis; Lcm, L. comus; Lcp, L. capensis; Lce, L. coreanus; Lci, L. corsicanus; Lcs, L. castroviejoi; Lfl, L. flavigularis; Lgr, L. granatensis; Lha, L. hainanus; Lin, L. insularis; Lmn, L. mandschuricus; Loi, L. oiostolus; Lot, L. othus; Lpg, L. peguensis; Lsn, L. sinensis; Lsx, L. saxatilis; Ltb, L. tibetanus; Ltl, L. tolai; Ltm, L. timidus; Ltw, L. towsendi; Nnt, Nesolagus netscheri; Ntm, N. timminsi; Ocn, Oryctolagus cuniculus; Opr, Ochotona princeps; Pfr, Pentalagus furnessi; Pmr, Poelagus marjorita; Pcr, Pronolagus crassicaudatus; Prn, Pronolagus randensis; Prp, Pronolagus rupestris; Psa, Pronolagus saundersiae; Rdi, Romerolagus diazi; San, Sylvilagus andinus; Saq, S. aquaticus; Sau, S. audubonii; Sbc, S. bachmani; Sbr, S. brasiliensis; Sdc, S. dicei; Sfl, S. floridanus; Sgb, S. gabbi; Smn, S. minensis; Snt, S. nuttallii; Sob, S. obscurus; Spl, S. palustris; Srb, S. robustus; Ssp, Sylvilagus sp.; Str, S. transitionalis; Tst, Tamiasciurus hudsonicus.