Introduction

Because of their conservative morphology, lagomorphs are notoriously difficult, to dissemble into biologically realistic evolutionary entities. ^{Bachman (1837}:282) notably stated that, “many of the species so greatly resemble each other in many particulars that the student in natural history has sometimes been greatly perplexed in deciding on the exact species referred to by authors.” Forty years later, ^{Allen (1877}), in listing examined specimens of Sylvilagus palustris- currently understood to be circumscribed to the southeastern United States (western limit to Mobile Bay)- listed specimens from Veracruz and Yucatán, México, as belonging in that species. It was not until the skull had been removed for examination that Allen instead admitted that the specimen from Veracruz should belong to its own discrete species, S. truei [= S. gabbi truei] (Allen, 1890b:192), noting that “the single record from so remote a point [i. e., from Florida] as Mirador, México, has of late seemed open to serious [taxonomic] question” thereby first remarking on congruence between geographic features and taxonomy of Sylvilagus. Taxonomic decisions at the time were routinely undertaken -with few but notable exceptions- based on external appearance. Since that time, increasingly detailed analyses have been undertaken, and an expanding tool chest of morphological characters have successively been employed to more accurately distinguish among lagomorph taxa (^{Baird 1857}; ^{Gray 1867}; ^{Lyon 1904}; ^{Nelson 1909}; ^{Thomas 1913}; ^{Hummelinck 1940}; ^{Hershkovitz 1950}; ^{Hall 1951}; ^{Palacios et al. 1980}; ^{Ruedas 1998}, 2017; Ruedas et al. 2017).

Philip ^{Hershkovitz’s 1950} treatise in particular, nominally focused on Colombia but in fact covering most, if not all, of South America, stands apart as the first attempt at a comprehensive treatment of the lagomorphs of any continent, albeit closely followed by E. R. ^{Hall’s 1951} synopsis of North American lagomorphs. The vast scope, both geographic and taxonomic, of Hershkovitz’s work meant that, years later, the taxonomy he proposed for Central and South American cottontails remained in force. For example, for Central American S. brasiliensis, ^{Hall (1951}, ¹⁹⁸¹) showed no changes relative to the scheme of Hershkovitz (1950). ^{Cabrera (1961}) similarly had few taxonomic changes in either S. “brasiliensis” or S. “floridanus” (both sensu lato) of South America, although S. nigronuchalis^{Hartert, 1894}, the oldest available name for South American taxa considered to be subsumed within S. floridanus, was inexplicably omitted from Cabrera’s treatment.

However, and notwithstanding its eminent worth, the passage of time has revealed that some errors made their way into ^{Hershkovitz’s 1950} treatise. As ^{Musser et al. (1998}:10) pointed out with particular respect to oryzomyine rodents, parts of some of Hershkovitz’s revisions could represent an “unfortunate example of taxonomic revision undocumented by specimens or other data and one that misleadingly simplified a complex reality”. In the case of the treatment of South American cottontails, for example, Hershkovitz did not examine all the pertinent holotypes, and when he did, it is unclear how carefully he scrutinized key morphological characters that could have led to a more accurate reflection of the underlying biological reality (^{Ruedas 2017}). In fact, Hershkovitz’s conclusion (1950:327) that his “review shows [S. brasiliensis and S. floridanus] to be the only recognizably valid species of leporids indigenous to South America” could not, in retrospect, have been further from the mark, given the recognized presence of a much larger number of species of Sylvilagus: at least 12 in the “brasiliensis” group alone (Ruedas et al. 2019).

In the present work, we began by questioning the taxonomy of individuals in the genus Sylvilagus from Costa Rica. Costa Rica, at 51,100 km², covers only 0.034 % of the land surface of the Earth, but with over 230 species present of terrestrial mammals, contains approximately 4 % of the World’s known mammal species: 121 times more than expected by strict proportionality between area and biodiversity. Insofar as cottontails (Sylvilagus) are concerned, that is reflected in the presence of three recognized taxa (^{Hall 1951}, 1981; ^{Mora 2000}; ^{Ruedas and Salazar-Bravo 2007}; ^{Rodríguez-Herrera et al. 2014}): S. g. gabbi (J. A. ^{Allen, 1877}), S. dicei^{Harris, 1932}, and S. floridanus costaricensis Harris, 1933. In describing S. f. costaricensis, ^{Harris (1933)} undertook comparisons of that taxon with S. f. aztecus (J. A. Allen, 1890) and S. f. hondurensis^{Goldman, 1932}. Goldman in turn, in his description of S. f. hondurensis, undertook comparisons between that taxon and S. f. chiapensis (^{Nelson, 1904}).

We accordingly undertook comparisons of taxa in the Sylvilagus floridanus group present in Costa Rica and the region in order to better ascertain their taxonomic identity. The Costa Rican-and indeed, Central American-taxa of Sylvilagus remain inadequately described, let alone diagnosed. We therefore undertook a detailed analysis of cranial and dental anatomy of Costa Rican taxa of S. floridanus within the broader context of their current nominal identification to species, by undertaking comparisons using all the pertinent holotypes: of S. floridanus (J. A. ^{Allen, 1890}): those of the species and subspecies listed above, and that of the geographically proximal S. f. yucatanicus (^{Miller, 1899}), thereby enabling us to robustly define the species of S. floridanus complex in Costa Rica and adjacent areas. Identification of species is, we believe, critical to generating phylogenetic trees that bear any semblance to the reality of life, because accurate trees can only result from the combination of adequate taxon sampling with sufficient data. Otherwise, one is left with what ^{Coddington and Scharff (1996}:139) so trenchantly remarked: “A fully resolved tree that makes no sense is still nonsensical.”

Materials and methods

Specimens. Specimens examined are listed in Appendix 1, with their original taxonomic designation as well as current taxonomy, localities (georeferenced insofar as possible), repository, and collection number. For geographic and taxonomic reasons, as described above, we chose to focus on the following taxa: Sylvilagus f. floridanus, S. f. costaricensis, S. f. hondurensis, S. f. aztecus, S. f. chiapensis, S. f. yucatanicus, S. gabbi, S. dicei, and S. brasiliensis surdaster (^{Thomas, 1901}). Sylvilagus b. surdaster was included because, although the type locality is in Ecuador (Esmeraldas Prov.; Río Bogotá, Carondelet; ca. 1° 07’ 27’’ N, 78° 45’ 45’’ W, ca. 20 m), and there would be scant probability of conspecificity, it is the most proximal lowland taxon affine to S. brasiliensis broadly writ and the name brasiliensis has previously been used for Costa Rican lowland rainforest rabbits following ^{Hall (1981}).

Morphological data: mensural characters. We measured 37 craniodental morphological variables. Terminology of cranial characters and features generally follows ^{Wible (2007}), and ^{Ruedas (1998}); measurements were defined by ^{White (1987}) and ^{Ruedas (1998}, ²⁰¹⁷), and were extensively detailed and illustrated in ^{Ruedas et al. (2017)}. Mensural characters included: GLS, greatest length of skull; POSTORB, width of postorbital constriction; BROSTR and DEPROSTR, breadth and depth (height) of rostrum; BBRAIN, breadth of braincase; ZYGO1, greatest width across the masseteric spine; ZYGO2, zygomatic breadth; LZYGO, length of zygomatic arch; NASALL, greatest length of nasal bone; NASALW, greatest width across left and right nasal bones; I2P2, least alveolar length of I2-P2 diastema; P2M3, greatest alveolar length of P2-M3 toothrow; HBRAIN, height of braincase; HBULLA, height of bulla; CONDL, condylopremaxillary length of cranium; LPALFOR, WPALFOR, length and width of incisive foramina; PALONG, palatal length; PALBRDGE, greatest anteroposterior dimension of palatal bridge; BASIOC, anteroposterior length of basioccipital; WIDBULL, width of auditory bulla; ANTBULL, anteroposterior length of auditory bulla, from the most anterior projection of the ectotympanic to the most posterior point between the occipital and the paracondylar processes of the exoccipital; INTBD, least breadth across the basioccipital between the ectotympanic bones; OCCOND, width across the occipital condyles; INTBOC, length between the posteriormost edge of the palatal bridge and the suture between the basioccipital and basisphenoid bones; CHOANA1, breadth of nasopharynx; CHOANA2, breadth of alisphenoid constriction; MASTOID, greatest breadth across the mastoid exposure of the petrosal; DEPZYGO, least anteroposterior length across the maxillary bone at the base of the masseteric spine on the maxillary portion of the zygomatic arch; IP3, least alveolar length of i-p3; MANDEP, depth of mandibular body; P3M3, greatest alveolar length of p3-m3; HMAND, height of the mandible; HPTT, distance from ventral aspect of angular process (labial to pterygoid shelf) to most dorsal aspect of pterygoid tuberosity; BCON, length of condyloid process; WCON, breadth of articular facet of condyloid process; LMAND, length of mandibular body.

Statistical analyses were carried out using the Statistical Analysis System (SAS) software, version 9.4 (2002-2012; ^{SAS Institute 1988a}, ^1988b), generally following ^{Ruedas (1995}, 1998); significance in all analyses was set at α = 0.05. Due to the paucity of specimens available, little could be made to determine presence or extent of sexual dimorphism in the taxa examined, although sexual dimorphism has been reported in measurements of Sylvilagus (^{Orr 1940}) and could affect results of multivariate analyses (^{Reyment et al. 1984}; ^{Marcus 1990}) given the small intraspecific sample sizes of the present study (Appendix 1). Univariate statistics (mean, standard deviation) were calculated using the UNIVARIATE procedure of SAS. Analysis of variance was carried out using the GLM procedure, enabling the MEANS routine with option REGWQ, which uses the Ryan-Eynot-Gabriel-Welsch multiple range test, and controls for Type I error (^{Day and Quinn 1989}). A principal component analysis (procedure PRINCOMP) was carried out on the covariance matrix of log-transformed normalized measurement values data. Such a posteriori grouping methods are preferred by us over a priori grouping methods (multiple range tests, discriminant analyses) because there is no prior hypothesis as to the putative identity of specimens examined. These data further are useful to examine ontogenetic growth patterns, which in the sample covariance matrix can be construed as the dispersion of points along the major long axis of each sample, with the first eigenvector representing Huxley’s allometric equation (^{Voss et al. 1990}). We used the broken stick method of ^{Frontier (1976}) as implemented by J^{ackson (1993}) to assess the significance of each principal component’s eigenvalue; broken-stick distributions for principal component eigenvalues were generated using the “broken.stick” function of R (v. 3.3.1; ^{R Core Team 2016}).

Dental characters. Drawings of p3 were made by tracing from photographs taken using a Canon EOS 30D digital camera mated to a Canon MP-E 65 mm f/2.8 1-5X Macro Photo lens, or a Canon EOS 6D mated to the same lens or an AmScope CA-CAN-SLR-III camera adapter for microscopes, shooting either through a camera tube on a binocular dissecting microscope or an ocular tube with the ocular removed, also on a binocular dissecting microscope, as made available by the collections housing the specimens under consideration. Among leporids, p3 generally constitutes the most informative dental element for taxonomic and systematic purposes (^{Dalquest 1979}; ^{Dalquest et al. 1989}; ^{Hibbard 1963}; ^{Palacios and López Martínez 1980}; ^{Ruedas 1998}; Ruedas et al. 2017; ^{White 1987}, ¹⁹⁹¹; ^{White and Morgan 1995}; ^{Winkler and Tomida 2011}). Discrete characters were deemed the most important in this particular research; accordingly, resulting figures were oriented and scaled to the same size in linear dimensions to carry out size-independent comparisons of interspecific characters. Characters considered follow the standard terminology of Palacios and López Martínez (1980), were described in Appendix I of ^{Ruedas (1998)} and illustrated here (Figure 1) with some modifications from Ruedas et al. (2017) in that all cusps are identified by incorporating features from ^{López Martínez (1974}, ¹⁹⁷⁷, ¹⁹⁸⁰, ¹⁹⁸⁹), ^{López-Martínez et al. (2007)}, and ^{Angelone and Sesé (2009}). Additional characters useful in distinguishing among lagomorph species were extracted from ^{Palacios (1996)} and ^{Palacios et al. (2008)}. The LSID for this publication is: urn:lsid:zoobank.org:pub:601C073B-6DFA-421E-8B4B-F7F44BF62D3F.

Figure 1 Standard nomenclature for dental features of Recent leporid lagomorphs’ third lower premolar (p3, top) and second upper premolar (P2, bottom), adapted from Figure 1 of Palacios and ^{López Martínez (1980}:62), and expanded from ^{Ruedas et al. (2017}) in identifying all cusps by incorporating features from López Martínez (1974, 1977, 1980, 1989) and ^{Angelone and Sesé (2009}). The term “anterior loph,” preferred herein, was used interchangeably with “trigonid” by ^{Hibbard (1963}). López-Martínez et al. (2007) considered only the caudal portion of the anterior loph of pm3 to constitute the trigonid, with the rostral portion (anterior lobe) instead collectively constituting the anteroconids. 

Results

Statistical analyses of morphology. Univariate statistics (means ± standard deviation, minimum-maximum) for the variables measured in each individual taxon (represented in certain taxa only by the holotypes or, in the case of S. gabbi, by the lectotype) are shown in Table 1. Also shown in Table 1 are the results of the Ryan-Einot-Gabriel-Welsch multiple range test. Thirty of the 37 characters examined showed some level of significance in discriminating among groups of individuals or taxa. This proportion (81.1 %) is markedly higher than the two characters that would be expected to differ significantly by chance alone with significance set at α = 0.05. However, some of the variables that are significantly different among taxa do not discriminate into distinct groups (e. g., depth of rostrum, mastoid breadth, length of mandibular toothrow, length of mandible). Similarly, most of the significantly different groups displayed a great deal of overlap. The one consistent result obtained from the analysis is that S. f. yucatanicus is immoderately larger than remaining taxa in almost all characters. That taxon differs significantly from all taxa but S. dicei in breadth of braincase, and from all other taxa in breadth of incisive foramina; it also has the longest skull of any Sylvilagus species examined for the present study, and beyond statistical significance (^{Moyé 2006}; ^{Wasserstein and Lazar 2016}; ^{Wassertstein et al. 2019}), does not overlap with the GLS of any of the remaining Sylvilagus taxa.

The results of the principal component analysis, carried out on the covariance matrix of a reduced set of natural log-transformed variables (n = 22; reduced as a compromise to embrace as many specimens as possible while maintaining as many measurements as possible), are shown in Figure 2 and Tables 2 and 3. The principal component analysis accounts for 15.4 % of the overall variance. Principal component 1 accounts for 36.2 % of that variation, with PC 2 accounting for 14.9 %; PCs 1-7 jointly account for >80 %, and 1-10 for >90 %. Just over half (50.5 %) of the variation in PC 1 is accounted for by only five of the 22 characters: width of bulla (13.5 %), width of incisive foramina (11.6 %), width of nasal bones (11.3 %), length of nasal bones (8.5 %), and breadth of rostrum (5.6 %). Remaining characters each contribute less than 5 % to the variation in his principal component.

Figure 2 shows the great deal of overlap among most taxa in the floridanus species group. Within the limits imposed by a reduced number of samples, the major axis of dispersion for points in these taxa is primarily along principal component 1, which in this instance is a size component. The major axis of dispersion has been shown to be associated with age-correlated growth (^{Voss et al. 1990}). In the particular instance of our analysis, this is borne out by the relative homogeneity of the magnitude of the eigenvector scores for PC 1 (Table 3): while some variables have eigenvector scores that are somewhat low (maxillary toothrow, interbullar distance, breadth of braincase) or somewhat high (width of bulla, width of nasal bones, width of incisive foramina), the remaining characters are fairly homogeneous, and eigenvectors average 0.200 ± 0.08. The homogeneity of the eigenvectors exhibited in PC1 is not evident in PC 2 through 6 or subsequent principal components (Table 3). For PC2, these average 0.041, but the standard deviation jumps to 0.209, with eigenvectors ranging from −0.456 (width of bulla) to 0.625 (breadth of nasopharynx). Subsequent principal components show similar trends with respect to standard deviation, maxima, and minima, of the characters’ eigenvector scores. Such lack of homogeneity in eigenvector scores usually is associated with shape-based, rather than size-based variation.

In Figure 2, principal component 2 (14.9 % of the total variation) distinguishes primarily between the floridanus group Sylvilagus species and remaining species, including S. gabbi and S. dicei. Only two characters contribute well over half (59.9 %) of the variation to this principal component: breadth of nasopharynx (39.1 %) and width of bulla (20.8 %).

Figure 2 Graphical results of the Principal Component Analysis undertaken on the correlation matrix of the reduced set (n = 22) of natural log-transformed variables. Letter codes as follows, a : S. f. aztecus; b : S. f. boylei; c : S. f. costaricensis; d : S. dicei; e : S. g. messorius; f : S. f. floridanus; g : S. g. gabbi; h : S. f. hondurensis; i : S. incitatus; m : S. f. margaritae; n : S. f. connectens; p : S. f. chiapensis; r : S. f. russatus; s : S. f. superciliaris; t : S. g. truei; u : S. f. daulensis; y : S. f. yucatanicus. Where the labelled polygon encloses its same designation letter (e. g., a , S. aztecus, or f , S. floridanus), the enclosed letter shows the location of the holotype in the first two dimensions of multivariate space; otherwise, letters refer to holotype (e. g., m , S. f. margaritae, or y , S. yucatanicus. 

The results of the principal component analysis reinforce the suggestion derived from the multiple range test that S. f. yucatanicus is exceptionally distinct from remaining taxa examined. That taxon is markedly separated in principal component 1 from remaining individuals examined (Figure 2), this despite the fact that we undertook natural log-transformation of the variables in order to minimize the effects of size. Width of incisive foramina is the second most important character in PC1, contributing to 11.6 % of the variation in that PC. Breadth of braincase in contrast only contributes to 1.1 % of the variation in PC 1.

Notwithstanding the informative nature of the exploratory principal components analysis, we acknowledge that said analysis is not without issues. Application of the broken stick method to assess the significance of eigenvalues suggested that only the first two principal components contained meaningful information. These two components cumulatively accounted for 51.1 % of the variation. Because the overall PCA accounted for 15.4 % of the variance, the result is that only 7.9 % of the morphological variation is accounted for in the PCA as implemented in the present study. It is possible that more judicious selection of variables may have influenced the analysis one way or another (e. g., selecting only those variables found to be significant in the multiple range test). We chose however to maintain the variables employed rather than cherry-pick the data. Our PCA results underscore that the morphological conservatism manifested in craniodental mensural variables throughout the genus Sylvilagus-and indeed, in other lagomorph genera-is not readily tractable to these morphometric analyses, although the analyses do have certain illuminative properties.

Taxonomic identity of Sylvilagus floridanus costaricensis ^{Harris, 1933}.

Analysis of morphological data. To ascertain the taxonomic identity of S. f. costaricensis, we undertook comparisons between this taxon and all other pertinent regional taxa, as noted in the introduction. Figures 3 4, 5 show the dorsal, ventral, and lateral perspectives, respectively, of the focal taxa: as one might expect from the results of the principal component analysis described above, observed differences among the various taxa are subtle. Nevertheless, they are present and telling. Notwithstanding, one obvious difference between the taxa under consideration is in greatest length of skull. In this character, our sample of S. f. floridanus show sexual dimorphism: in males, the mean in mm ± SD (min-max) is 70.8 ± 1.1 (69.4 to 72.1), whilst in females, it is 72.7 ± 0.9 (71.3 to 73.8); t ₉ = 2.9125, P < 0.0172, δ _means = 2.0 mm, 95 % CI = 0.4-3.5 mm. However, our sample of adults of S. f. hondurensis includes only one female (AMNH 126205); remaining individuals are either unknown (AMNH 123378) or males (n = 7). Our comparisons in measurements are therefore made grouping the sexes. Between S. f. floridanus and S. f. hondurensis, the respective data are 72.0 ± 1.4 (69.4 to 73.8), versus 76.4 ± 1.0 (74.6 to 77.6), t ₁₈ = 7.6512, P < 0.0001, δ _means = 4.4 mm, 95 % CI = 3.1-5.6 mm. The holotype of S. f. costaricensis, at 76.3 mm, is congruent with the mean of S. f. hondurensis. The same pattern obtains, albeit without sexual dimorphism in S. f. floridanus (P = 0.1142), in breadth of skull at the zygomatic spine: 33.9 ± 0.6 (33.0 to 35.0), versus 35.5 ± 0.4 (34.9 to 36.1), t ₁₅ = 6.033, P < 0.0001, δ _means = 1.7 mm, 95% CI = 1.1-2.2 mm.

The region of the frontonasal suture, and the shape of the latter, is a character that has been used extensively in previous taxonomic studies of lagomorphs. For a selected subset of the specimens employed herein, that feature is shown in Figure 6. The specimens in the top row all are S. f. floridanus collected contemporaneously; these all show the posterodorsal process of the premaxilla extending caudad of the frontonasal suture (even with the terminus in USNM 76711), along with a short, marked intrusion of the frontal bone extending between the posterodorsal process of the premaxilla and the nasal bone. The caudally projecting posterodorsal process of the premaxilla is apparent in Central American taxa only in S. g. gabbi. The nasal bones themselves are significantly smaller in S. f. floridanus than in S. f. hondurensis: means in mm ± SD (min-max) are respectively 30.5 ± 1.0 (28.4 to 31.9), versus 34.8 mm ± 0.9 (33.2 to 36.3), t ₁₇ = 9.2366, P < 0.0001, δ _means = 4.3 mm, 95 % CI = 3.3-5.2 mm. The holotype of S. f. costaricensis is congruent with S. f. hondurensis in nasal bone length (35.3 mm), and in morphology in that the posterodorsal process of the maxilla is retracted rostrally relative to the caudal terminus of the nasal bone. One might expect that because of the longer GLS, the Central American taxon would naturally have a longer nasal bone. However, the Pearson product-moment correlation coefficients between GLS and NASAL suggest that this is not necessarily the case: for S. f. hondurensis, R = 0.611 (R ² = 0.373, P = 0.108), whereas for S. f. floridanus R = 0.753 (R ² = 0.567, P = 0.007); the holotype of S. f. costaricensis is almost identical in these two measurements to S. f. hondurensis AMNH 126203. We therefore predict that given larger sample sizes, S. f. costaricensis will be more closely allied to the pattern displayed by S. f. hondurensis.

A corollary of the shorter nasal bone in S. f. floridanus is that that bone does not extend as close to the orbit in S. f. floridanus as in S. f. costaricensis and S. f. hondurensis (Figure 3). Measured from the most posterolateral point of the nasal, the distance to the caudalmost point in the notch between the antorbital process and the frontal bone is 6.3 mm in S. f. floridanus, 3.5 mm in S. f. hondurensis, and 3.8 mm in S. f. costaricensis.

Figure 3 Dorsal views of the crania of the Central and South American taxa under consideration herein (current nomenclature), scaled to the same greatest length of skull. A: Sylvilagus f. floridanus, holotype, AMNH 1890/1155 ( ♀ ), greatest length of skull (GLS): 71.9 mm; B: S. f. costaricensis, holotype, UMMZ 65232 ( ♀ ), GLS: 76.3 mm; C: S. f. hondurensis, holotype, USNM 257062 ( ♂ ), GLS: 74.7 mm; D: S. f. aztecus, holotype, AMNH 3116/2438 ( ♂ ), GLS: 72.9 mm; E: S. f. chiapensis, holotype, USNM 75953 ( ♀ ), 78.9 mm; F: S. f. yucatanicus, holotype, USNM 37772 ( ♀ ), GLS: 81.1 mm; G: S g. gabbi, lectotype, USNM 11371/37794 ( ♂ ), GLS: 70.9 mm; H: S. dicei, holotype, UMMZ 64043 ( ♀ ), 77.3 mm; I: S. brasiliensis surdaster, holotype, MNH 1901.6.5.16 ( ♀ ), GLS: 72.7 mm. 

In the holotype of S. f. costaricensis, there is a small intrusion of frontal bone, the nasopremaxillary process of the frontal, separating the caudal tip of posterodorsal process of the premaxilla from the caudal tip of the nasal bone (measured from the tip of the posterodorsal process of the premaxilla, right: 3.3 mm, left: 4.5 mm). This intrusion is absent from the holotype of S. f. hondurensis and largely absent from examined specimens in this taxon, although some (e. g., AMNH 123378, Figure 6) have a minute manifestation of this feature. The frontonasal suture also may vary in shape, being either parallel with a transverse plane starting laterally then angling rostrally to meet the opposite nasal bone at the medial plane, or on an approximate diagonal plane in a caudo-lateral to rostromedial direction. Sylvilagus f. floridanus displays the former, whereas S. f. costaricensis represents the latter condition; in this character, S. f. hondurensis is more similar to S. f. floridanus.

Other characters of the dorsal aspect are somewhat more shrouded. Pitting in the parietal and frontal bones has for example been employed as a character in distinguishing between taxa (^{Wible 2007}; ^{Ruedas et al. 2017}; ^{Ruedas 2017}). However, there is a thin layer of tissue covering this portion of the skull of the holotype of S. f. costaricensis that, despite its slenderness, obscures this character. Similarly, the angle of the suture between the parietal and supraoccipital is somewhat descending ventrally from external to medial direction in S. f. floridanus, but is horizontal or ascending in S. f. hondurensis. However, it is not clearly visible in the holotype of S. f. costaricensis.

Figure 4 Ventral views of the crania of the Central American taxa under consideration herein, scaled to the same width. Specimens shown are the same as in Figure 3 and disposed in the same order. 

From a lateral perspective (Figure 5), the length of the zygomatic arches of S. f. floridanus differ significantly with little overlap in size from those of S. f. hondurensis: 31.0 ± 0.8 (29.3 to 32.2) v. 33.0 ± 1.1 (31.0 to 34.3), t ₁₈ = 4.621, P = 0.0002, δ _means = 2.0 mm, 95 % CI = 1.1-2.9 mm. The zygomatic arch of S. f. costaricensis, at 32.2 mm, is at the upper limit of those of S. f. floridanus, but is firmly ensconced within those of S. f. hondurensis. The relative brevity of the zygomatic arch of S. f. floridanus gives it a more robust dorsoventral appearance than those of S. f. costaricensis and S. f. hondurensis; however, vertical depth of the zygomatic arch does not differ significantly among the taxa: 5.3 ± 0.5 (4.3 to 6.0) v. 5.0 ± 0.3 (4.6 to 5.3), t ₁₇ = 1.3771, P = 0.1863, δ _means = 0.3 mm, 95 % CI = −0.21-0.7 mm; the holotype of S. f. costaricensis has a zygomatic depth of 5.3 mm, congruent with either taxon. As in length of nasal bones, this likely is a manifestation of the differences in GLS, given that S. f. floridanus and S. f. hondurensis have almost identical zygomatic length relative to GLS: 43.0 % and 43.2 %; 42.1 % in S. f. costaricensis. Also as in the nasal bones, however, the length of the zygomatic arch is significantly correlated with GLS in S. f. floridanus (R = 0.761, R ² = 0.579, P = 0.006), but not in S. f. hondurensis (R = 0.384, R ² = 0.148, P = 0.307).

Figure 5 Lateral views of the crania of the Central American taxa under consideration herein, scaled to the same width. Specimens shown are the same as in Figure 3 and disposed in the same order. Inverted for consistency are: S. f. costaricensis and S. f. aztecus. The latter also was not taken on a completely lateral plane, making the profile appear more dorsoventrally bowed than it is in reality. 

Analysis of dental morphology. Substantial and substantive differences are exhibited in the dental morphology the taxa under consideration herein (Figure 7). In the tooth most commonly used to discriminate among species of lagomorphs, lower premolar 3, S. f. costaricensis differs from S. f. floridanus in several key features: the anteroflexid is relatively deep and U-shaped, with a narrow constriction on the rostral surface, whereas in S. f. floridanus, the anteroflexid is broadly open and V-shaped; S. f. hondurensis displays a condition similar to S. f. costaricensis. Other Central and South American comparator taxa examined here display a more complex pattern on the rostral surface of pm3, with multiple anteroflexids or, if single, with a complex internal structure (e. g., S. f. chiapensis). In S. dicei, the rostral architecture of pm3 is of such complexity that a lingual anteroconid is identifiable as a region only, rather than as a distinct feature of the tooth.

The central angle, an almost universal feature of the lagomorph rostral hypoflexid, is present as a singular inflection in S. floridanus, but as an unusual double inflection in S. f. costaricensis. The central angle is indistinct in S. f. chiapensis because of the complexity of the enamel pattern, and possibly double in the lectotype of S. g. gabbi. Sylvilagus brasiliensis surdaster displays a very weak central angle. The caudal surface of the hypoflexid is relatively smooth (labial portion) to somewhat crenulate (lingual portion) in S. f. floridanus versus highly complex and strongly crenulate in S. f. costaricensis. Other regional taxa display the range from similarly crenulate morphologies (S. dicei, S. f. yucatanicus) to somewhat less crenulate (S. f. aztecus, S. f. chiapensis), to completely smooth (S. b. surdaster).

While S. f. floridanus definitively does not exhibit a paraflexid (being instead convex), there is a slight inflection in that portion of pm3 of S. f. costaricensis. Sylvilagus f. hondurensis has a concavity at the base of the anteroconid that we likewise interpret as a paraflexid, as does S. b. surdaster. Otherwise, this surface of the tooth is relatively featureless from slightly convex (S. f. yucatanicus) to slightly concave (S. f. aztecus, S. f. chiapensis).

Figure 6 Detail of the left frontonasal suture and posterodorsal process of premaxillary bone in selected individuals, all scaled to the same anteroposterior length for consistency. A: Sylvilagus f. floridanus, AMNH 1890/1155 ( ♀ , holotype; greatest length of nasal, in straight line from facialmost to caudalmost point: 30.4 mm); B: USNM 70870 ( ♀ ; 30.6 mm); C: USNM 76711 ( ♂; 29.4 mm); D: USNM 77113 (subadult ♂ ; 22.4 mm; note the difference in proportions of dimensions); E: USNM 77114 ( ♀ ; 31.1 mm); F: USNM 77115 ( ♂; 31.4 mm); G: S. f. hondurensis, USNM 257062 ( ♂ , holotype; 34.2 mm); H: S. f. hondurensis, AMNH 126146 ( ♂ ; 35.0 mm); I: S. f. costaricensis, UMMZ 65232 ( ♀ , holotype; 35.3 mm); J: S. f. yucatanicus, USNM 37772 ( ♀ , holotype; 37.2 mm); K: S. g. gabbi, USNM 11371/37794 ( ♂ , lectotype; 27.9 mm); L: S. dicei, UMMZ 64043 ( ♀ , holotype; 33.6 mm). Key features include: caudal terminus of the posterodorsal process of premaxilla relative to the caudal terminus of the nasal bone, and absence, presence, and rostral extent of process on frontal bone extending between posterodorsal process of premaxilla and posterolateral margin of nasal bone (nasopremaxillary process of frontal bone). 

In PM2 of S. f. floridanus, the hypoflexus is marked by a slight depression, barely demarcating mesial from distal hypercones. In S. f. costaricensis, there is a distinct, deep, U-shaped hypoflexus in PM2, resulting in distinct mesial and distal hypercones. The area of PM2 between postcone and poststyle PM2 in S. f. floridanus is convex, with no trace of a metaflexus. In contrast, S. f. costaricensis has a small but distinct inflection marking the metaflexus.

The first upper incisor, although generally neglected as featureless among lagomorphs, also is distinct between the two taxa: in S. f. costaricensis, lingual and labial cusps are subequal in height relative to the rostral groove demarcating them; S. f. hondurensis is almost identical in the morphology of its I1. In contrast, the lingual cusp of I1 in S. f. floridanus is distinctly expanded rostrally relative to the labial cusp.

Figure 7 Crown views of the enamel structure of lower premolar 3 (upper two rows), upper premolar 2 (middle two rows), and first upper incisor (lower two rows) for the Central American taxa under consideration herein. Specimens in each triplet are, A: Sylvilagus f. floridanus, holotype, AMNH 1890/1155 ( ♀ ); B: S. f. costaricensis, holotype, UMMZ 65232 ( ♀ ); C: S. f. hondurensis, holotype, USNM 257062 ( ♂ ); D: S. f. aztecus, holotype, AMNH 3116/2438 ( ♂ ); E: S. f. chiapensis, holotype, USNM 75953 ( ♀ ); F: S. f. yucatanicus, holotype, USNM 37772 ( ♀ ); G: S g. gabbi, lectotype, USNM 11371/37794 ( ♂ ); H: S. dicei, holotype, UMMZ 64043 ( ♀ ); I: S. brasiliensis surdaster, holotype, MNH 1901.6.5.16 ( ♀ ). Some images were rotated horizontally in order for all perspectives to be the same; all images are scaled to the same width so as to show differences in proportion rather than in size. In each image, rostral is to the top of the figure, labial is to the right of the figure. 

An additional, and unusual, feature is present in I2 of S. f. costaricensis. This tooth is invariably small, cylindrical, with a circular cross section in every species of Sylvilagus we have examined to date. However, in S. f. costaricensis, I2 is roughly triangular in cross section, with the base caudal and apex rostral, and has two distinct grooves on the caudal aspect of the tooth (Figure 8). The only other taxon of Sylvilagus that we have examined for this study to display these characters is S. f. hondurensis.

Taxonomic conclusion: identity of Sylvilagus floridanus costaricensis. In light of the foregoing analyses, particularly those based on cranial and dental characters, it is clear that the differences between S. f. costaricensis and S. f. floridanus are interspecific in nature insofar as taxa of Sylvilagus are concerned. As described above, the skulls differ significantly in magnitude in a number of measurements; they also differ significantly in a number of cranial and dental characters. However, S. f. costaricensis are not distinct from S. f. hondurensis in the same characters. Most significantly, both taxa share two unique synapomorphies: a triangular cross section to I2, which is marked by two grooves on its caudal facies. We therefore consider that S. f. costaricensis are not distinct from S. f. hondurensis at the species level. Sylvilagus floridanus hondurensis was described by Edward A. ^{Goldman on 30 July 1932}; S. f. costaricensis by William P. ^{Harris on 28 June 1933}. As a consequence, the name hondurensis has priority. Until a greater number of specimens are available for examination of population level and broader extent of geographic variation, there are sufficient differences between the two taxa-for example, the comparative extent and degree of crenelation of the caudal aspect of the pm3 hypoflexid-that we recommend the prudent cause of action to keep both names, as Sylvilagus hondurensis hondurensis E. A. Goldman, 1932, and S. hondurensis costaricensis Harris, 1933.

Figure 8 Occlusal perspective of the first and second right upper incisors of S. f. costaricensis (left) and S. f. hondurensis (right). Arrows mark the two grooves on the caudal aspect of I2. Note the unusual triangular cross section of I2, rather than the almost universal condition for Sylvilagus of a circular cross section for this tooth. 

Discussion

We consider our study foundational to any future regional or focused taxonomic study of biogeography, evolution, and phylogeny of cottontails. Revolutions in the practice of taxonomy and phylogenetics have led to a more nuanced understanding of species delimitation and, as a result, of species boundaries. ^{Ruedas et al. (2017}) noted that there is a lack of cohesion between philosophical and operational approaches to species; as in that work, we apply what ^{Sangster (2014}) called “methodological introgression” of species concepts applied in an operationally coherent manner to “discover, describe, and order into our classification system” (^{Mayden 1997}:387) the individuals within, or constituting, the species category, independent of the properties of the species category. We used previously (^{Ruedas et al. 2017}) an integrative approach to species delimitation (sensu ^{Padial et al. 2010}; ^{Schlick-Steiner et al. 2010}) as implemented by ^{Naomi (2011}). This approach, using a morphological character set vastly expanded over that of ^{Hershkovitz (1950}), resulted in hypotheses of taxonomic species in Sylvilagus that reflected the underlying biological reality imposed by abiotic criteria such as elevation, temperature, and precipitation regimes, soils, etc., as well as the effects of those abiotic factors on vegetation, which ultimately is reflected by the species inhabiting the ecosystems under consideration. While there have been controversies regarding the application of, for example, the phylogenetic species concept to particular instances (e. g., ^{Groves and Grubb 2011} vs. ^{Zachos et al. 2013}; ^{Zachos 2015}), the integrative approach yields coherent and biologically relevant taxonomic hypotheses: a single widespread species of Sylvilagus (S. “brasiliensis” sensu Linnaeus 1758) distributed from the Atlantic to the Pacific coasts of South America, from 0 to >5,000 m in elevation, and from Veracruz, México, in the north, to Argentina in the south is neither coherent, nor biologically realistic. The taxonomic hypotheses we propose herein for S. floridanus follow from Allen’s hypothesis that geography, while not the ultimate arbiter of taxonomy, nevertheless strongly affects species limits: “Hence the single record from so remote a point […] has of late seemed open to serious question” (^{Allen 1890}:192). The biogeographic breaks in Central and South America, reflected in the taxonomy of numerous taxa, are likewise reflected in Sylvilagus. In South America, rivers have been implicated in speciation events in small mammals (^{da Silva and Patton 1998}; ^{Matocq et al. 2000}; ^{Patton et al. 2000}), primates (^{Wallace 1852}; ^{Boubli et al. 2015}), and birds (^{Naka and Brumfield 2018}) alike. Sylvilagus are similarly affected by vicariant effects. In the instance of Sylvilagus, the effects of strong ecological change brought about by the xeric conditions at the Isthmus of Tehuantepec also appear important.

From a biogeographic perspective, the patterns of speciation revealed by our taxonomic framework are congruent with those of other taxa. For example, Bassariscus astutus is restricted to the north and west of the Isthmus of Tehuantepec, and its sister species B. sumichrasti, while somewhat overlapping the range of B. astutus in coastal Guerrero and Oaxaca, México, largely is restricted to the east and south of the isthmus. Similarly, taxa in the Reithrodontomys sumichrasti species complex (Rodentia: Cricetidae: Neotominae) show an analogous distribution and hypothesized relationships (^{Hardy et al. 2013}) as we propose here for Sylvilagus. In the case of mice of the genus Habromys (Rodentia: Cricetidae: Neotominae), six of the seven species in the genus are restricted to the north and west of the isthmus, and only one species, H. lophurus, is restricted to the south and east of the isthmus (^{León-Paniagua et al. 2007}). This pattern of sister taxa of mammals exclusively distributed to one or the other side of the Isthmus of Tehuantepec is a repeating evolutionary and biogeographic motif (^{Sullivan et al. 2000}; ^{Rogers et al. 2007}).

One result of the integration of distinct data streams to assess taxonomic relationships is the stark difference in taxonomic information content that is brought about by using morphometric (continuously variable measurement data) versus discrete character data. Our principal components analysis (Figure 2) shows that there is substantial overlap in morphology among the distinct taxa of Sylvilagus when these are subjected to morphometric analysis. Of note, the principal components analysis is an a posteriori test, thus there is no prior hypothesis imposed on the ensuing result. In contrast, an a priori test such as a discriminant function analysis essentially “forces” the output to conform to the a priori hypothesis (i. e., predict group-species-membership) because it describes a function that will distinguish among the predefined samples groups (i. e., presumptive taxa). As a result, a posteriori tests are preferable in taxonomy because they do not impose a hypothesis on the data, rather the results derived from the data are a reflection of the presumptively true nature of the underlying taxonomic reality. In the present instance, however, the two statistically significant principal components only accounted for 7.9 % of the morphological variation among the groups. That is to say, conversely, that 92.1 % of the mensural variation went unaccounted for. Thus, either a posteriori or a priori tests would be on tenuous grounds in terms of establishing-or even testing-a robust taxonomic hypothesis, no doubt because of the morphologically conservative, or strongly homomorphic nature of cranial morphology in Sylvilagus, and indeed, in Leporidae in general. Because of this, and based on the results of our morphometric analysis, taxa clearly distinguished in the analysis (e. g., gabbi, dicei, yucatanicus) are hypothesized to be definitively distinct; however, taxa in our sample that overlap in multivariate space are not definitively demonstrated to be the same, i. e., subject to Type II error. It is in these circumstances that inspection of character data becomes increasingly valuable: assessment of discretely variable characters in morphologically conservative taxa, particularly when such characters may be discretely distinct in morphometrically indistinguishable groups, can result not only in identification and discrimination of different taxa but also in the possibility of inferring evolutionary relationships among the groups or taxa in question. Character data (qualitative) can be useful for identifying and classifying organisms, while morphometric data (quantitative) may under certain circumstances be useful for identifying organisms, as well as for studying the physical (mensural) characteristics of organisms and their variation. Excessive reliance on either, particularly morphometric data, may result in erroneous taxonomic hypotheses.

Unanswered, however, remains the question of: why are there so many species of Sylvilagus present in Costa Rica? We hypothesize that the present biodiversity is a combination of the ecological heterogeneity of Costa Rica, along with its location. We have previously documented, using molecular approaches (^{Ruedas et al. 2017}), that there were multiple invasions of South America by Sylvilagus. Some of the remaining biodiversity of Costa Rican Sylvilagus may be essentially remnants of these multiple invasions: taxa that resulted from populations that remained in place as other populations continued to expand the dispersal front. As remnant populations, their conservation therefore becomes ever more imperative.

Taxonomic Conclusions. On the basis of the foregoing, we recognize the following taxa in Central America south of the Isthmus of Tehuantepec to the Panama-Colombia border: Sylvilagus dicei, S. gabbi, and S. hondurensis.

Sylvilagus hondurensis

Honduras cottontail

Sylvilagus floridanus hondurensis^{Goldman, 1932}:122. Type locality, “From Monte Redondo, about 30 miles northwest of Tegucigalpa, Honduras (altitude about 5,100 feet [1554 m]).” The village of Monte Redondo lies at ca. 860 m rather than, as indicated by Goldman, at 1,554 m. Roads lead NW from Monte Redondo to higher elevations. The 1,554 m contour on a road emanating from Monte Redondo is at ca. 14° 18’ 42’’ N, 87° 18’ 24’’ W. We speculate that Goldman referred to the higher elevations today contained within the Reserva de Vida Silvestre Corralitos (Francisco Morazán, Honduras), just NW from the village of Monte Redondo. Holotype: USNM 257062.

Sylvilagus floridanus costaricensis^{Harris, 1933}:3. Type locality, “from Hacienda Santa Maria, Province of Guanacaste, Costa Rica, altitude 3,200 feet” (975 m). The Hacienda Santa María ranger station, inside Guanacaste National Park is located at 10° 45’ 52’’ N, 85° 18’ 11’’ W, 844 m, thus corresponding fairly closely with Harris’ description. Holotype: UMMZ 65232.

Sylvilagus yucatanicus

Yucatan cottontail

Lepus aquaticus: ^{Allen, 1877}:365 (part). Not Lepus aquaticus^{Bachman, 1837}. Allen noted that “In the collection are quite a number of specimens from the provinces of Vera Cruz and Yucatan in Southern México. These differ from specimens from Mississippi and Louisiana in no very marked degree.” He later revised his opinion (^{Allen 1890b}) and transferred these specimens to Lepus sylvaticus [= S. floridanus].

Lepus sylvaticus aztecus: ^{Allen, 1890}:191, from “Merida, Yucatan”; not Allen 1890:188, from “Tehuantepec City”.

Lepus floridanus yucatanicus^{Miller, 1899}:384. Type locality, “Merida, Yucatan” (correctly spelled “Mérida, Yucatán” by ^{Hall 1951}:159). Holotype, USNM 11441/37772.

Sylvilagus floridanus yucatanicus: ^{Lyon, 1904}:336. Name combination.