Introduction

Understanding the relationships among evolutionary lineages is critical to estimating species diversity at varying spatial scales, reconstructing the evolutionary history of taxa, delineating ecological communities, and in making informed conservation decisions (^{Crozier 1992}; ^{Faith 1992}; ^{Crandall et al. 2000}; ^{Sinclair et al. 2005}; ^{Chave et al. 2007}). Molecular data are increasingly used for evaluating relationships among species, identifying potential species-level clades, and identifying so-called cryptic species and thus can have significant impact on our understanding of evolutionary relationships (^{Sinclair et al. 2005}; ^{Beheregaray and Caccone 2007}; ^{Bickford et al. 2007}; ^{Mort et al. 2015}).

With the continual improvement of molecular techniques and analyses, and broader sampling of natural populations, our understanding of phylogenetic relationships is often in flux. A group that has proven particularly difficult to delineate with traditional morphological characters are the species of spiny pocket mice of the genus Heteromys (Rodentia: Heteromyidae; see ^{Anderson 2015} and references therein). ^{Goldman (1911)}, in the first revision of the genus, recognized 13 species of Heteromys dividing them into two subgenera: Heteromys containing 12 species and Xylomys with a single species. This author further recognized the subgenus Heteromys as comprising two distinct species groups with the H. desmarestianus group containing eight species, including the first named Heteromys and most widely distributed species, H. desmarestianus (^{Gray, 1868}). Recent research based on mitochondrial DNA indicates, however, that the lowland dry forest spiny pocket mice that were long recognized as a separate and sister genus, Liomys, are paraphyletic with respect to the species of Heteromys, thus should either be recognized as species of Heteromys or as another generic level clade (see ^{Anderson et al. 2006}; ^{Hafner et al. 2007}; ^{Anderson and Gutiérrez 2009}; and references therein). The current trend is to consider all species as belonging in the genus Heteromys.

Heteromys desmarestianus has remained a recognized taxon through several revisions (^{Hall 1981}; ^{Rogers and Schmidly 1982}; ^{Williams et al. 1993}; ^{Patton 2005}; ^{Rogers and González 2010}), although new species have since been recognized (^{Anderson and Jarrín-V 2002}; ^{Anderson 2003}; ^{Anderson and Timm 2006}; ^{Anderson and Jansa 2007}; ^{Anderson and Gutiérrez 2009}). Recently, ^{Rogers and González (2010)} suggested four additional clades within H. desmarestianus should be recognized. This research focuses on one of those four proposed clades that is located within the Caribbean lowlands of Costa Rica.

Heteromys desmarestianus, as currently defined, is common and widespread, ranging from southern Mexico to Colombia (^{Reid 2009}). This species is found in evergreen and semideciduous forests, from sea level to high elevation cloud forests (^{Timm et al. 1989}; ^{Reid 2009}). In Central America’s Caribbean lowlands, the forest spiny pocket mouse is difficult to study because populations are often found at low densities (^{Fleming 1974}; ^{Timm et al. 1989}), and anthropogenic disturbances often have negative impacts on density and species diversity (Romero, pers. obs.). Based on molecular evidence from mitochondrial and nuclear DNA of three individuals from Caribbean lowlands of Costa Rica, ^{Rogers and González (2010)} suggested that these individuals may actually represent a separate species from what is recognized as H. desmarestianus. Herein, we test across multiple sites in the lowlands, if individuals from the Caribbean lowlands of Costa Rica are genetically distinct from what is recognized as H. desmarestianus and how this population(s) and others of the H. desmarestianus species complex are related to each other. In order to build a better understanding of the species diversity in this lineage and to test the hypothesis that there is greater diversity than is currently recognized we, herein, evaluate the relationship within the lineage currently recognized as the species H. desmarestianus.

Materials and methods

We trapped mice in several locations throughout the Caribbean lowlands of Costa Rica from 2007-2010 (Figure 1, Table 1). The Caribbean lowlands have similar ambient temperature throughout, but annual precipitation can vary from 2,400 to 4,800 mm per year (^{McClearn et al. 2016}). Our localities (Table 1) ranged in elevation and size of forested area; our individuals from the highest elevation were from the Berlin property (Destierro) ranging from 210 to 280 masl. Samples from Berlin also represented our southern-most sample. Our northern-most samples were from the Refugio Nacional de Vida Silvestre Mixto Maquenque, close to the Costa Rica-Nicaragua border on the Río San Juan (Figure 1).

Figure 1 Map of localities for all specimens used in the study. Localities 29-38 represent specimens from Costa Rica’s Caribbean lowlands. Symbols correspond to lineages depicted in the maximum likelihood tree (Figure 2). H. anomalus from Venezuela not depicted herein. Specific data regarding localities can be obtained from Appendix 1. 

We used Sherman live traps (8 cm × 9 cm × 23 cm; H. B. Sherman Traps, Inc., Tallahassee, FL) placed at ground level and baited with cracked corn, oats, and mixed bird seed. Traps were checked daily, and when an individual was caught, a toe was removed with surgical scissors and immediately placed in 95 % ethanol. All vials with tissue and ethanol were stored frozen within hours of collection. Voucher specimens of both complete specimens and toe samples are deposited at the University of Kansas Natural History Museum, Lawrence, Kansas. This project was undertaken with the approval of the University of Kansas Institutional Animal Care and Use Committee. All animal handling protocols were in accordance with the guidelines of the American Society of Mammalogists (^{Sikes et al. 2016}).

Laboratory procedures—116 samples from 10 sites in the Caribbean lowlands were used (Table 1) for genetic comparisons. Tissues were soaked in deionized water for one hr prior to beginning the digestion process. Standard digestion and DNA extraction were conducted following the protocol for mouse tails in ^{Sambrook et al. (1989)}. The mitochondrial cytochrome-b (cytb) gene was amplified in full using the primers 765 and 766 (^{Bickham et al. 2004}). Polymerase chain reaction (PCR) was performed using 50 μl reactions of the following reagents: 5 μl of 10× buffer, 5 μl of 10× MgCl₂, 5 μl of 10× solution of dNTP, 0.5 μl of Taq DNA polymerase, 5 μl of a 10× solution of each primer, 25 μl of deionized water, and 1-2 μl of extracted DNA. Thermal cycle conditions consisted of initial heating at 94°C for 3 min, then 36 cycles of denaturation at 94°C for 30 s, annealing at 50°C for 30 s, and extension at 70°C for 2.5 min. PCR products were purified using the QiAquick PCR purification kit (Qiagen, Valencia, CA) and were subsequently used in standard sequencing reactions using Big Dye version 3.0 (Applied Biosystems, Foster City, CA). Sequencing reactions were cleaned using Sephadex spin columns and analyzed with an ABI 3100 automated genetic analyzer (Applied Biosystems, Forster City, CA). Sequence data were manually aligned using Sequencher v. 4.9 (Gene Codes Corporation, Ann Arbor, MI). We used the complete cytb gene, and all flanking regions were discarded prior to phylogenetic analysis.

To expand our dataset, we used cytb sequences of H. desmarestianus available from GenBank (^{Benson et al. 2013}). We incorporated 74 individuals representing samples from throughout the range of the species, including specimens from near the type locality of H. desmarestianus, Cobán, Guatemala. Heteromys anomalus, H. australis, and H. nelsoni were used as outgroups (Appendix 1). We aligned all sequences with Muscle v.3.8.31 (^{Edgar 2004}) implemented in Jalview 2.8 (^{Waterhouse et al. 2009}).

Phylogenetic analysis—Phylogenetic relationships were obtained by performing a maximum likelihood (ML) analysis. We estimated models of molecular evolution using jModelTest v.2.1.1 with the corrected Akaike information criterion test (^{Guindon and Gascuel 2003}; ^{Darriba et al. 2012}). We used GARLI v. 2.0 (^{Zwickl 2006}) for ML analyses, using two independent search runs, with a maximum of five million generations each. Support values were calculated using bootstrap with 500 replications in GARLI, and results visualized and edited in FigTree v.1.4 (^{Rambaut 2007}).

Results

The aligned data set comprises 1,142 characters of which 738 were constant, 335 characters were parsimony-informative, and 69 variable characters were parsimony-uninformative. The model of DNA substitution inferred from jModeltest 2.1.1 is TIM2+I+G.

The ML tree topology (Figure 2) shows two highly supported lineages for all individuals currently considered H. desmarestianus. One clade comprises all samples from the Caribbean lowlands of Costa Rica and had very strong (99 %) bootstrap support. The other lineage comprises all of the H. desmarestianus sequences obtained from Belize, Costa Rica (in part), El Salvador, Guatemala, Honduras, México, and Nicaragua and also have strong bootstrap support (83%). Within the clade containing samples from the Caribbean lowlands of Costa Rica (not shown in Figure 2), the individuals from two sites, Berlin and Maquenque, formed clades supported by strong bootstrap support (90 % and 86 %, respectively). Two individuals from the eight sampled at Maquenque were placed elsewhere within the tree, and thus the nonexclusive nature of the branching pattern complicates lower level population patterns from these data. Maquenque is biologically quite interesting being in the floodplain of the Río San Juan and our continued studies there are elucidating other unexpected patterns with other rodent species.

Figure 2 Maximum likelihood tree of our Heteromys dataset, based on cytb sequences. The maximum likelihood tree is collapsed for visual clarity, and shows the results with bootstrap values. These results indicate that individuals from Costa Rica’s Caribbean lowlands harbor unique mitochondrial lineages that could reflect reproductive isolation. Symbols in tree are used in Figure 1 to show the geographic range of lineages. Tree is rooted with H. golmani. 

Specimens collected near Baja Verapas, Guatemala (GU646966, GU646967, GU646968; Appendix) ~30 km from the type locality, Cobán, cluster with specimens from Belize, El Salvador, Honduras, and parts of México, but not with specimens from the Caribbean lowlands of Costa Rica (Figure 2). The specimen that clustered closest to Cobán from our Costa Rican lowland dataset was from the Costa Rica-Nicaragua border region ~ 850 km away. Genbank accession numbers of new sequences are reported in the Appendix.

Discussion

The results of this study strongly support the hypothesis that what is currently called Heteromys desmarestianus in Central America and northern-most South America comprises two very distinct clades, one being found in the Caribbean lowlands of Costa Rica and the other comprising all other populations. Our results show a clear geographic pattern; individuals currently considered H. desmarestianus in the Caribbean lowlands of Costa Rica harbor distinct mtDNA haplotypes from individuals considered H. desmarestianus found elsewhere in the Neotropics, including other areas in Costa Rica (Figure 2). The specimens from near the type locality of H. desmarestianus, Cobán, Guatemala form a well-supported clade with specimens from southern México, Belize, El Salvador, and Honduras. Specimens from western Costa Rica and one specimen from western Nicaragua all form another distinct clade.

Costa Rica is only ~51,000 km², yet its variable topography and climate result in diverse habitats with unique flora and fauna (^{Janzen 1983}). Currently, four main mountain ranges divide the country into the Pacific and Caribbean sides. These mountain ranges span southeast to northwest, and are of diverse ages and origins (^{Anderson and Timm 2006}). Extending from western Panama to northern Costa Rica, the Cordillera de Talamanca, Cordillera Central, and Cordillera de Tilarán form an expansive mountain range with peak elevations of over 3,000, 2,500, and 2,000 masl, respectively. The Cordillera de Guanacaste is the northernmost range in Costa Rica, and is comprised of several isolated volcanoes, with passes of ~500 to 700 masl in elevation that connect the Pacific and Caribbean sides (^{Anderson and Timm 2006}). The historical and current topography of these mountain ranges probably shaped the diversification and speciation patterns in the flora and fauna observable today.

The family Heteromyidae originated on the North America continent (^{Wood 1935}; ^{Schmidly et al. 1993}), and fossil remains for the subfamily Heteromyinae are known from the Pliocene, Pleistocene, and Holocene (^{Rogers 1990}). ^{Rogers (1990)} estimated that the major groups within this subfamily diverged ~12 to 13 mya, yet the historical events that produced the H. desmarestianus group are largely unknown. The geologic history of the Caribbean, and Central and South America has been a debated topic (^{Bartoli et al. 2005}; ^{Montes et al. 2012a}, ^2012b, ²⁰¹⁵; ^{Bacon et al. 2015}; ^{O’Dea et al. 2016}), but it is thought that islands of volcanic origin between Central and South America may have allowed faunal exchanges prior to the formation of a permanent land bridge (^{Bartoli et al. 2005}; ^{Woodburne 2010}; ^{O’Dea et al. 2016}). The time of the emergence of a permanent Panamanian land bridge is disputed, and estimates range from 2 to 7 mya (^{Montes et al. 2012b}; but see ^{Bacon et al. 2015}). Because of the widespread distribution pattern of the H. desmarestianus group, a hypothesis similar to the one suggested for other rodent groups has been proposed for this clade (^{Patterson and Pasqual 1972}; ^{Baskin 1978}; ^{Simpson 1980}; ^{Rogers 1990}; ^{Almendra and Rogers 2012}; ^{Pine et al. 2012}); It is thought that considerable radiation occurred in the Miocene and Pliocene throughout Central America, with a subsequent entry to South America via the Panamanian land bridge (^{Rogers 1990}; ^{Schmidly et al. 1993}).

Our results herein document a broad distribution of this distinct and unrecognized Costa Rican lowland lineage. The northern-most samples within this clade are from the Costa Rica-Nicaraguan border, while the southern-most are ~86 km southeast of there. Unfortunately, little is known about spiny pocket mice in the lowlands of Nicaragua and southern Costa Rica, and we are unable to demarcate northern and southern boundaries of this mtDNA lineage. In terms of elevation, the lowland specimens came from forests that ranged in elevation from ~22 to 280 masl. In our analysis, we included a single GenBank sequence from Cerro Honduras in Parque Nacional Braulio Carrillo. The park, along with privately owned reserves and biological stations, is part of a continuously forested transect that expands from the lowlands at La Selva Biological Station and reaches elevations > 2,700 masl. Although we do not have specific data on the elevation from which this particular specimen came, our results do show that this higher elevation specimen is a member of the clade with H. desmarestianus proper from throughout Central America and does not group with our samples from the lowlands, including specimens from the nearby La Selva Biological Station. While our results suggest two distinct clades, we are unable to delineate at this time limits of their specific elevational range, or if there are areas of overlap or hybridization that await discovery.

Our results expand upon, compliment, and confirm ^{Rogers and González (2010)}, who used both cytb and nuclear data, and identified three individuals from the Caribbean lowlands of Costa Rica as a potential candidate species. Although useful to characterize species that are difficult to establish based on morphological data, DNA sequence data do have limitations, particularly when a single marker is used (^{Farias et al. 2001}; ^{Rogers and González 2010}). Now that we provide more extensive sampling of individuals from the Caribbean lowlands of Costa Rica, we recommend that future studies determine if population structure based on nuclear markers correlates with the distinctive mtDNA lineage of the Costa Rican lowland. Further investigation focused on nuclear DNA is also important because mitochondrial DNA and nuclear DNA can be discordant (^{Lack et al. 2010}, ^{Bernardo et al. 2019}). This can result in distinct mitochondrial DNA lineages within a population or species that are not supported by nuclear DNA.

The diversity of rodent communities in the Caribbean lowlands of Costa Rica have been vastly understudied and we believe underestimated, in part, because of low densities resulting in low trap success (Romero, pers. obs.). Consequently, the lack of data and specimens has hindered our understanding of the basic phylogenetic relationships and biogeographic patterns of species in the area. Other widespread rodent species have been found to hold similar patterns reported herein, where individuals from the Caribbean lowlands of Costa Rica are genetically distinct and potentially new species (Timm, unpublished data). These data suggest that there may be significant cryptic diversity in the lowlands, and that more phylogenetic studies should include samples from this region to identify potential biogeographic patterns for rodents in the Neotropics. This information is necessary not only to understand phylogenetic relationships, but also to have a grasp on the patterns and levels of diversity for the area, and make large-scale conservation decisions based on this information. We believe that our results, in conjunction with future studies that aim to identify and delineate diversity in the H. desmarestianus species complex, and the relationships between these species, will allow for a greater understanding of the historical events leading to speciation in this group.

Clearly much remains to be learned about the diversity of these widespread and common rodents that are considered keystone species in the Neotropics.