Introduction

Remote oceanic islands and archipelagos are biologically simpler than continental regions and therefore provide ideal geographical and historical settings for the study of colonization, adaptation, speciation, and diversification of species. Islands have long been recognized as natural models for the study of evolutionary processes (^{Parent et al. 2008}; ^{Losos 2010}; ^{Rodrigues and Diniz-Filho 2016}; ^{Román-Palacios and Wiens 2018}). Compared to continental regions, islands are more ideal places to observe and interpret patterns of evolution due to their geographic isolation, small size, fewer numbers of species, and a high degree of endemism. The development of ecological and evolutionary systems can be directly observed on volcanic islands as we are now better able to date the timing of their emergence above the ocean surface as blank slates for colonization and the timing of subsequent evolutionary diversification (^{Losos and Ricklefs 2009}; ^{Hendriks et al. 2019}).

The Galápagos Islands are a young oceanic and volcanic archipelago resulting from the eastward passage of the Nazca plate over a hotspot, at a rate of 59 km/My, located in the Pacific Ocean at approximately 960 km west of the coast of South America. It is composed of 13 major islands larger than 10 km², six smaller islands, over 40 islets with official names and many smaller unnamed islets and rocks, for a total of approximately 8,000 km² of land spread over 45,000 km² of water (^{Snell et al. 1996}; ^{Parent et al. 2008}, ^{Geist et al. 2014}, ^{Harpp et al. 2014}). The age of the islands increases moving eastward with the oldest islands located towards the southeast of the archipelago. The present islands date from up to 3.5 to 4 million years ago (Ma) for the eastern islands of Española and San Cristóbal, respectively, to 60,000 years ago for Fernandina Island (^{Geist et al. 2014}).

Despite their tropical climate, the Galápagos Islands have been the stage of surprisingly few animal diversifications compared with other Pacific tropical island groups. Among vertebrates, the absence of amphibians and the virtual absence of mammals are particularly striking and nearly unique among terrestrial island ecosystems (^{Parent et al. 2008}). However, ^{Román-Palacios and Wiens (2018)} showed that the Galápagos archipelago drove faster rates of speciation and diversification in tanagers and tortoises, at least three times higher than in other related lineages of birds and tortoises inhabiting different islands.

The rodents of the Galápagos Islands are the only terrestrial mammals that have naturally colonized the islands and diversified within the archipelago. They belong to the Neotropical rice rat family Cricetidae, subfamily Sigmodontinae, tribe Oryzomyini. Two genera of rodents are currently known to be present on the islands: Aegialomys (^{Weksler et al. 2006}) and Nesoryzomys (^{Heller 1904}). There are two commonly recognized and described species in the genus Aegialomys: A. xanthaeolus the type species from the mainland and A. galapagoensis, formerly known from San Cristobal Island but has not been collected since their initial capture by Darwin in 1835 and is presumed extinct there (^{Patton and Hafner 1983}; ^{Dowler et al. 2000}). In addition, the subspecies A. galapagoensis bauri, extant in Santa Fé Island, was previously considered a separate species but has been most often considered synonymous with A. galapagoensis (^{Cabrera 1961}; ^{Musser and Carleton 1993}, ²⁰⁰⁵; ^{Weksler et al. 2006}; ^{Weksler and Percequillo 2011}; ^{do Prado and Percequillo 2018}). The genus Nesoryzomys comprises two extinct and three extant species. The two extinct species are N. indefessus from Santa Cruz and Baltra Islands and N. darwini from Santa Cruz Island. The extant species include N. swarthi from Santiago Island, and N. narboroughi and N. fernandinae from Fernandina Island (^{Dowler and Carroll 1996}; ^{Dowler et al. 2000}). N. narboroughi has sometimes been synonymized with N. indefessus (^{Heller 1904}; ^{Musser and Carleton 1993}, ²⁰⁰⁵), but in this manuscript we retain both N. indefessus and N. narboroughi as different species as recommended by ^{Dowler (2015)}. The giant rice rat, Megaoryzomys curioi, is known only from subfossil remains from Santa Cruz Island and is not known to have a mainland representative (^{Patton and Hafner 1983}). It is possible that its extinction occurred prior to human settlement of the archipelago.

To date, the only systematic study to include nearly all of the species in the Galápagos, both extinct and extant, is ^{Patton and Hafner (1983)}. They did not analyze N. fernandinae (^{Hutterer and Hirsch 1980}), which was described as a new species after their research was in press (^{Dowler et al. 2000}). Based on a variety of data sets, including morphology, anatomy, protein electrophoresis and chromosome number and morphology, they suggested that: a) Nesoryzomys should be recognized at the generic level, nevertheless its origin is ambiguous, b) there were at least two independent colonizations of the islands, with Nesoryzomys representing an early arrival at 3 to 3.5 Ma, followed considerably later by Aegialomys as late as a few hundred to a thousand years ago, c) both Aegialomys taxa from the islands derived from A. xanthaeolus of the coastal Peruvian river valleys, d) N. narboroughi, N. swarthi and N. indefessus should be considered as races of a single species; and e) A. galapagoensis and A. bauri should be considered conspecific.

^{Dowler et al. (2000)} performed one of the most recent field surveys during which they discovered a population of N. fernandinae; this allowed the first opportunity to describe the appearance of this species, which was previously known only from skeletal remains (^{Hutterer and Hirsch 1980}) from Fernandina Island. They also found a viable population of N. swarthi on Santiago Island, which had previously been presumed extinct. These specimens represent the first endemic rodents taken on Santiago Island since the type series was collected in 1906 (^{Orr 1938}), and a single partial skull was found in 1965 (^{Peterson 1966}). However, recent molecular phylogenetic studies included only A. xanthaeolus, N. narboroughi and N. swarthi. ^{Weksler (2003}, ²⁰⁰⁶⁾ and ^{Weksler et al. (2006)}, using a nuclear exon and morphology, found that N. narboroughi and N. swarthi are monophyletic, and they are closely related to A. xanthaeolus. ^{Pine et al. (2012)}, ^{Leite et al. (2014)}, and ^{Machado et al. (2014)}, using morphology and mitochondrial and nuclear genes, supported the results of ^{Weksler (2003}, ²⁰⁰⁶⁾ and the results of ^{Weksler et al. (2006)}.

Outside of what was learned from these studies, little is known about the endemic rodents of the Galápagos Islands; however, it is clear that they are critically threatened by invasive species and human activities. The goal of the present study was to undertake the first study including comprehensive sampling of all four extant endemic rodent species (A. galapagoensis, N. narboroughi, N. swarthi and N. fernandinae) inhabiting the archipelago, to provide a dated phylogeny, and to elucidate the population genetics status of each species. This information will help to elucidate the evolutionary history of these island taxa and synthesize information about evolution and biogeography at scales that span remote islands, archipelagoes and continents.

Methods

Sample collection. We obtained tissue samples (liver, kidney) from museum specimens deposited at Angelo State Natural History Collections (ASNHC) at Angelo State University and the Museum of Vertebrate Zoology at the University of California, Berkeley (MVZ). Additional samples were from ear biopsies from animals released at the collection site. We sampled 159 individuals of the Aegialomys and Nesoryzomys genera inhabiting the Galápagos Islands, A. galapagoensis (n = 43), N. narboroughi (n = 49), N. swarthi (n = 43), and N. fernandinae (n = 24; Figure 1, Appendix 1), and two samples of A. xanthaeolus from Ecuador to elucidate its relationship with A. galapagoensis. We used Pseudoryzomys simplex and Oligoryzomys microtis as outgroups for the phylogenetic analyses (sequences downloaded from GenBank, accession numbers AY863422.1 and AY863420.1, respectively).

Figure 1 Sampling localities of A. galapagoensis (red dots; n = 43), N. narboroughi (pink dots; n = 49), N. swarthi (blue dots; n = 43) and N. fernandinae (green dots; n = 24) in the Galápagos Islands. 

DNA isolation and mitochondrial gene amplification. We extracted DNA from tissues using the DNeasy® Blood and Tissue Kit (QIAGEN, Inc., Valencia, CA) by cutting approximately 20 μg of tissue into several small pieces and following the manufacturer’s protocol. We amplified the mitochondrial control region (D-loop) by polymerase chain reaction (PCR) using the primers designed for Oligoryzomys spp. (^{González-Ittig et al. 2002}). In rodents, the D-loop has been useful in phylogenetic analyses due to its elevated mutation rate, lack of recombination and maternal inheritance (^{Robins et al. 2014}). Also, due to its high mutation rate, this marker has been used to detect signatures of population structure at a scale of just a few kilometers (^{Hirota et al. 2004}; ^{Urgoiti et al. 2018}). The PCR reactions contained 14.85 μL DEPC H2O, 3.0 μL 10X Reaction buffer, 3.0 μL deoxynucleotide triphosphates (dNTPs, 2 mM of each), 1.5 μL each forward and reverse primer (10 μM), 3.0 μL of 0.1 % bovine serum albumin (BSA), 0.15 μL AmpliTaq DNA polymerase (Applied Biosystems, Foster City, CA), and 3 μL of DNA for a final reaction volume of 30 μL. We used 1 % agarose gels stained with ethidium bromide to visualize DNA extractions and to amplify products. Amplification parameters were as follows: initial step of 95 ºC (4 min), followed by 34 cycles of denaturation at 95 ºC (40 sec), annealing at 50 ºC (30 sec), extension at 72 ºC (90 sec), and a final extension at 72 ºC (10 min). Reaction products were purified using AMPure Magnetic Beads (Agencourt Bioscience, Beverly, MA). The Oligoryzomys D-loop primers were used in an initial sequencing run. More specific internal forward and reverse primers were designed using the sequence fragments obtained from this initial run. The internal primers GIF (5’ - CCACTACCAGCACCCAAAGCTG - 3’) and GIREV (5’- GGTTGTGTTGATTAATGATCC - 3’) were used in all sequencing reactions. Five microliters of cleaned PCR product were added to 4μL of ABI Prism Big Dye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems, Foster City, CA) and 1.0 μL of 1.6 μM internal primer in each sequencing reaction (GIF and GIFREV). Sequencing reaction conditions were: 96 ºC (1 min), 45 cycles of 96 ºC (30 sec), 58 ºC for GIF or 52 ºC for GIREV (15 sec), 60 ºC (4 min), followed by a final holding step of 4 ºC. All sequencing reactions were performed using GeneAmp® PCR System 9700 (Applied Biosystems, Foster City, CA). Final sequencing products were purified using Sephadex G-50 powder then dried in a vacuum centrifuge and stored at - 20 ºC. Sequences were re-hydrated with the addition of 5 μL of HiDi Formamide with 0.1 mM EDTA, denatured at 95 ºC (3 min) and sequenced with capillary action electrophoresis using SCE 9610 Genetic Analysis System (SpectruMedix, State College, PA).

Phylogenetic analyses and divergence times estimation. We cleaned and edited sequences using Geneious® 11.1.4 (https://www.geneious.com), and performed multiple sequence alignment using ClustalW v.2.1 (^{Larkin et al. 2007}) implemented in Geneious. The best evolutionary model of nucleotide substitution was estimated in jModelTest 2.1.1 (^{Guindon and Gascuel 2003}; ^{Darriba et al. 2012}) using the Akaike information criterion. A Bayesian Inference (BI) analysis was performed in MrBayes 3.2.6 (^{Huelsenbeck and Ronquist 2001}; ^{Ronquist and Huelsenbeck 2003}), run for 20 million generations sampling every 1,000 generations. Output parameters were visualized using Tracer v1.7.1 (^{Rambaut et al. 2018}) to check for convergence between runs, and the first 25 % of the trees were discarded as burn-in.

We used BEAST v2.5.2 (^{Bouckaert et al. 2014}) to estimate molecular dates of divergences under an uncorrelated lognormal relaxed molecular clock model. The time to the most recent common ancestor for the main lineages was obtained using Bayesian Markov chain Monte Carlo (MCMC) searches. We sampled trees and divergence dates for all nodes every 10,000 iterations for 50,000,000 generations. These analyses implemented the Yule speciation processes model and the randomly generated starting tree as priors. We used three calibration points. The first calibration was based on a biogeographical event: the origin of the Galápagos archipelago at 5 Ma (^{Geist et al. 2014}). ^{Machado et al. (2014)} found that the lineage leading to the endemic genus Nesoryzomys derived from an ancestor shared with the clade composed of Melanomys, Sigmodontomys and Aegialomys and split around 1.49 (95 % HPD : 0.26 to 3.23) Ma. We used the split between Nesoryzomys and Aegialomys as a second calibration point. The third calibration point was based on the split of the lineage leading to Pseudoryzomys around 2.58 (95 % HPD : 0.43 to 5.38) Ma (^{Machado et al. 2014}). We checked convergence statistics for effective sample sizes using Tracer v1.7.1 (^{Rambaut et al. 2018}). We used TreeAnnotator v2.5.2 (available in the BEAST package) to get a consensus tree with node height distribution after elimination of 25 % of trees as burn-in. We visualized MrBayes and Beast results using FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/).

We performed BI and Beast analyses including all samples per species (trees not shown). We chose to perform the phylogenetic and dating analyses using representative samples of each species to reduce the saturation effects. Populations-level analyses included all the samples.

Population analyses. We conducted population genetic analyses separately for each clade determined with the BI analysis, with the exception of A. xanthaeolus because of the small sample size. The number of haplotypes (H), nucleotide diversity (π) (^{Nei 1987}), haplotypic diversity (h), number of polymorphic segregating sites (S), singletons (S₁), parsimony informative sites (PIS), and the average number of nucleotide differences (K) were estimated using DNAsp v6.12 (^{Rozas et al. 2017}). We performed Fu’s Fs (^{Fu 1997}) and Tajima’s D (^{Tajima 1989a}) neutrality tests to evaluate whether data departed from a neutral model of evolution due to factors such as population bottleneck or sudden expansion. Statistical significance was determined using the coalescent simulator in DNAsp v6.12 (^{Rozas et al. 2017}) with 1000 simulations.

We used the distribution of the number of pairwise mutational differences among individuals, or mismatch distribution, to explore demographic patterns of populations using DNAsp v6.12 (^{Rozas et al. 2017}); graphical representation was made by means of the growth-decline model. Raggedness (r) index (^{Harpending 1994}) and R ₂ statistics of ^{Ramos-Onsins and Rozas (2002)} were calculated to analyze goodness of fit of a population expansion model using 1000 simulations in the same program. Populations at demographic equilibrium or in decline should provide a multimodal distribution of pairwise differences, whereas populations that have experienced a sudden demographic expansion should display a star-shaped phylogeny and a unimodal distribution (^{Tajima 1989b}; ^{Slatkin and Hudson 1991}; ^{Rogers and Harpending 1992}; ^{Harpending and Rogers 2000}). However, recent changes in population size may not be detectable in mismatch distribution analyses due to threshold effects, time lags, or earlier demographic events that may mask the effects of recent events (^{Rogers and Harpending 1992}; ^{Harpending and Rogers 2000}).

We analyzed the magnitude of historical demographic events by constructing Bayesian Skyline Plots (BSP) using BEAST v2.5.2 (^{Bouckaert et al. 2014}). This analysis infers population fluctuations over time by estimating the posterior distribution of the effective population size at specific intervals along a phylogeny (^{Drummond and Rambaut 2007}). Genealogies and model parameters were sampled every 10,000 iterations along 50,000,000 generations under a relaxed molecular clock, with 25 % of burn-in. Convergence statistics for effective sample sizes and demographic plots were visualized using Tracer v1.7.1 (^{Rambaut et al. 2018}). In comparison with simple parametric and older coalescent demographic methods, the smoother estimates and sensitivity of this method, together with credibility intervals, provide a realistic population size function and enable retrieval of more details than just summary statistics (^{Deli et al. 2016}).

To further investigate the genetic relationships of the haplotypes at the intraspecific level we constructed haplotype networks using the Median-Joining algorithm (^{Bandelt et al. 1999}) implemented in PopART v1.7 (^{Leigh and Bryant 2015}).

Results

Phylogenetic lineages and divergence time. We amplified an average of 643 bp of the D-loop gene from 159 individuals of the Aegialomys and Nesoryzomys genera inhabiting the Galápagos Islands (A. galapagoensis n = 43, N. narboroughi n = 49, N. swarthi n = 43, and N. fernandinae n = 24, Figure 1, Appendix 1) and two samples of A. xathaeolus from the mainland.

The TVM with gamma distribution ( + G) model was recognized as the best fitting model with the following parameters: base frequencies A = 0.3612, C = 0.2445, G = 0.1109, T = 0.2834; nst = 6; and rates = gamma with shape parameter (α) = 0.6020. The BI and Beast analyses including all the samples per species (trees not shown) confirmed the monophyly of each of the species. However, the relationship with the outgroups was not well resolved and the posterior probabilities were lower. This could be due to the high mutation rates that are inherent to the D - loop region which can result in genetic saturation. Distantly related taxa are often affected by saturation effects. When sequences in a multiple alignment have undergone multiple substitutions, the apparent distances largely underestimate the real genetic distances and the alignment is said to be saturated (^{Philippe et al. 2011}). In phylogenetics, saturation effects result in long branch attraction, decrease of phylogenetic information, and underestimation of observed divergence times (^{Wilke et al. 2009}, ^{Philippe et al. 2011}). We thus focused our subsequent results and conclusions on the analyses that included only representative samples of each species. The BI analyses showed a topology in which four main clades for the Galápagos species were recognized with high levels of support and one clade including the continental species (Figure 2). All of the Galápagos rodent species analyzed were monophyletic. Clade 1 corresponds to all individuals recognized as A. xanthaeolus, which is a sister group of Clade 2, which includes all representatives of A. galapagoensis. Clade 3 corresponds to all samples identified as N. narboroughi, which is the sister to Clade 4 and Clade 5, corresponding to N. swarthi, and N. fernandinae, respectively.

Figure 2 Phylogenetic tree of the four extant endemic rodents of the Galápagos Islands based on Bayesian Inference analysis of mtDNA D-loop sequence data. Numbers at nodes indicate support values of posterior probabilities. 

The analysis estimating the time to the most recent common ancestor (TMRCA) showed that the oldest divergence event corresponds at the split between the Nesoryzomys and Aegialomys genera, dated around 3.84 (95 % HPD : 2.91 - 4.88) Ma (Figure 3). The first split between a continental versus an island species, A. xanthaeolus and A. galapagoensis, occurred around 1.11 (95 % HPD : 0.37 to 2.11) Ma. Speciation within the genus Nesoryzomys started around 2.23 (95 % HPD : 1.32 to 3.12) Ma with the split between N. fernandinae and N. swarthi versus N. narboroughi, followed by the split between the two lineages composed of N. fernandinae and N. swarthi which occurred around 1.58 (95 % HPD : 0.91 to 2.42) Ma. The estimated dates of divergence for the main nodes and their highest posterior density values are shown in Table 1. According to these results, the main speciation events occurred since the early-middle Pliocene up to the Pleistocene. However, the diversification within each species started at the end of the Pleistocene around 525,500 years ago.

Figure 3 Dated phylogeny of the extant endemic rodents of the Galápagos Islands reconstructed from mtDNA D-loop haplotypes inferred from BEAST. The horizontal bars show the 95 % confidence intervals. Time-scale in millions of years ago. Dates and letters at nodes depict values calculated in Table 1. 

Demographic reconstruction. Genetic diversity and neutrality test per species are shown in Table 2. All of the endemic species of the Galápagos Islands showed high genetic diversity (Hd > 0.965). N. fernandinae has the highest number of unique haplotypes in proportion with the number of samples, and N. swarthi the lowest.

The applied neutrality test revealed significant deviations from mutation-drift equilibrium for all the species inhabiting Galápagos according to Fu´s F _s values, and only for A. galapagoensis using Tajima´s D. The negative values suggest recent population expansion events in these species (Table 2).

The statistical analyses of mismatch distribution showed unimodal distributions for A. galapagoensis, N. narboroughi and N. fernandinae, which also suggests a recent demographic expansion (^{Slatkin and Hudson 1991}, ^{Rogers and Harpending 1992}) or spatial expansion (^{Ray et al. 2003}, ^{Excoffier 2004}). Statistical analysis of the mismatch distribution r index and R ₂ were significant (Figure 4).

Figure 4 Mismatch distribution of pairwise differences of haplotypes for each of the extant species inhabiting the Galápagos Islands. Shown are observed (red lines) and expected (dark blue square dot lines) frequencies obtained under a model allowing for population size change. Raggedness (r) index and R ₂ statistics values are shown. Significance is indicated with a star (*). 

The Bayesian skyline demographic reconstructions showed a pattern of constant size of populations through time, followed by a recent and small increase (Figure 5). A. galapagoensis and N. fernandinae show a small growth period, which started around 10,000 and 13,000 years ago, respectively. N. narboroughi and N. swarthi show a more constant trend of population increase, starting around 50,000 and 90,000 years ago, with a rapid increase starting around 6,000 and 25,000 years ago, respectively. Skyline reconstruction showed that A. galapagoensis, N. narboroughi, N. swarthi and N. fernandinae have similar populations sizes.

Figure 5 Skyline plots for each extant species of rodents inhabiting the Galápagos Islands. Plots show posterior median (darker lines) and 95 % Bayesian credible intervals (lighter lines on the outside) of the effective population size. The Y axis is in logarithm scale. 

Overall, the network analyses including all haplotypes for each species show a few abundant haplotypes, with frequencies between two and six, and numerous unique ones for all the species (Figure 6). N. fernandinae showed more unique haplotypes, with two being the highest frequency observed for a haplotype. The four species show long branches with haplotypes that are highly differentiated from the other haplotypes. N. narboroughi and N. fernandinae, for which samples were collected from three and two localities in Fernandina Island, respectively, do not show that the distribution of the haplotypes follows any structure or differentiation. Only N. narboroughi has two shared haplotypes among localities. Though A. galapagoensis and N. narboroughi do not show networks with a star-like shape, both networks show haplotypes with many connections, suggesting recent populations or with recent demographic expansion. In contrast, the networks of N. swarthi and N. fernandinae may suggest older or more stable populations.

Figure 6 Median-Joining network for each species of the extant rodent of the Galápagos Islands, using mtDNA D-loop haplotypes. Circles sizes are proportional to the haplotype’s frequencies in the population. The lines along the branches connecting the haplotypes show the number of substitutions. 

Discussion

This is the first genetic study to include all extant endemic rodent species of the Galápagos Islands. We implemented different genetic analyses in order to elucidate the phylogenetic relationships among these rodents inhabiting the archipelago, as well as demographic history and relationships of their populations. This information is extremely important for conservation of these endemic species, given that there is very little known about their biology and ecology, and that they are considered vulnerable, mainly due to human activities and the introduction of invasive species. Furthermore, we contributed phylogenetic and demographic information as well as divergence estimates in order to form hypotheses regarding colonization of the islands and compare them with previous hypotheses.

Our phylogenetic study corroborated the monophyly of the genera Aegialomys and Nesoryzomys as ^{Patton and Hafner (1983)} and ^{Weksler (2003}, ²⁰⁰⁶⁾, proposed. However, those authors did not include all the extant species. We also corroborated the monophyly of the continental species A. xanthaeolus and the island species A. galapagoensis, N. naboroughi, N. swarthi and N. fernandinae. Our calibration results suggest that the main speciation events started during the Pliocene with the split between the genus Aegialomys and Nesoryzomys (3.84 Ma, 95 % HPD:2.9 to 4.88), which agrees with the time proposed by ^{Patton and Hafner (1983)}, using Nei’s methods, dated around 3 to 3.5 Ma; it differs from that proposed by ^{Machado et al. (2014)}, dated around 1.49 Ma (95 % HDP:0.26 to 3.23), and ^{Parada et al. (2013)} around 2.4 Ma (95 % HDP not available) during the Pleistocene. The difference between Machado et al.’s estimate and ours could be due to the genes used; they used IRBP and Cyt b, while we used D - loop. Also, despite the fact that our study and Machado et al.’s both used the origin of the Galápagos as a calibration point, we used the new date proposed by ^{Geist et al. (2014)} of 5 Ma, while Machado et al. used 4 Ma (^{Geist 1984}). We also included A. galapagoensis from the islands, while they only used A. xanthaeolus. Including the two species of the genus Aegialomys, and specifying an older date for the origin of the archipelago, resulted in an older date for the split between these genera, which is well supported with the paleogeographic and biogeographical information. Despite this discrepancy, there is substantial overlap in the credibility intervals of these two assessments.

The second major speciation event occurred within the genus Nesoryzomys with the split between N. narboroughi versus N. swarthi and N. fernandinae dated at 2.23 Ma (95 % HDP:1.32 to 3.12) during the Pleistocene. ^{Leite et al. (2014)} dated the split between N. narboroughi and N. swarthi around 2 Ma (95 % HDP not available) during the Pleistocene, supporting our results. The third major event occurred at 1.58 Ma (95 % HDP:0.91 to 2.42) during the Pleistocene with the split between N. swarthi and N. fernandinae. Finally, the speciation between A. xanthaeolus and A. galapagoensis occurred at 1.11 Ma (95 % HDP: 0.37 to 2.10), also during the Pleistocene.

Rice rats are the only terrestrial mammals that naturally colonized the Galápagos and diversified within the archipelago (^{Clark 1984}, ^{Parent et al. 2008}). There are some prior hypotheses about how these two genera colonized islands. ^{Patton and Hafner (1983)} suggested that Nesoryzomys is an old immigrant to the islands while the Aegialomys species is quite recent. They based their results on the high degree of morphological and biochemical (allozyme) distinctiveness of Nesoryzomys relative to other oryzomyines, including Aegialomys. This supported a more ancient origin and a single immigration to the islands, whether of pre-Nesoryzomys or Nesoryzomys form with subsequent radiation within the islands. In contrast A. galapagoensis is nearly morphological identical to the mainland A. xanthaeolus, which suggests an introduction to the islands within the last few hundred to thousand years, perhaps via pre-Columbian humans coursing the west coast of Perú (^{Patton 1984}; ^{Pine et al. 2012}). ^{Parent et al. (2008)} suggested that colonization events on the Galápagos occurred over the last 3 to 4 Ma during the existence of the present islands, and might have happened earlier when now sunken islands were above sea level. The presence of several drowned seamounts on the Carnegie Ridge east of the Galápagos (^{Christie et al. 1992}) suggests that earlier volcanic islands may have served as stepping stones for colonization by some of the terrestrial fauna. ^{Parent et al. (2008)} also mentioned that rice rat diversity is the result of one colonization event from the South American continent, where their close relatives inhabit, as it has been shown in other species of tortoises and lizards (i. e.,^{Caccone et al. 1999}; ^{Kizirian et al. 2004}; ^{Benavides et al. 2007}).

Our results suggest two colonization events to the Galápagos by the extant species. The first arrival by the ancestor of Nesoryzomys was probably a species of Sigmodontomys or Melanomys that originated from lower montane and lowland forest habitats in South America (^{Harris and Macdonald 2007}). This event should be dated from the early-late Pliocene boundary onwards, as proposed by ^{Leite et al. (2014)}. This was followed by a second arrival of Aegialomys from coastal Perú by rafting over recent historical times to the middle Pleistocene, giving rise to the Aegialomys of the Galápagos (^{Patton 1975}; ^{Steadman 1985}; ^{Hutterer and Oromí 1993}; ^{Weksler 2003}). It is well known that these species, and the most closely related ones, are excellent dispersers across salt water (^{Pine et al. 2012}), making colonization easier for them.

The fauna and flora of the Galápagos Islands is principally derived from western South America, Central America, and the Caribbean (^{Merlen 2014}). So, we assume that colonization from there to the archipelago during the Pliocene was possible because there is strong evidence (^{Christie et al. 1992}; ^{Werner et al. 1999}) that islands have been forming over the hotspot for at least twice as long as the age of the oldest islands, and perhaps as long as 20 Ma. These ancient islands are now seamounts east of the present-day Galápagos along the Carnegie Ridge on the Nazca Plate (^{Christie et al. 1992}) and northeast of the archipelago on the Cocos Plate on the Cocos Ridge (^{Werner et al. 1999}). There has been a “conveyor belt” of islands produced over millions of years, providing the potential for organic colonization during that time (^{Merlen 2014}).

Once the ancestors of Nesoryzomys and Aegialomys arrived to the Galápagos Islands, dispersal, colonization, speciation, and diversification were possible during the Pleistocene, in part, because the integrated area of the Galápagos Islands was much greater than today, land bridges existed between a few of the major islands, and many more minor islands and islets were exposed. In fact, there is a hypothesis that between 1 and 5 Ma, at least nineteen major Galápagos Islands existed but are currently submerged; these are in addition to the thirteen that exist today (^{Geist et al. 2014}). Later island fragmentation has led to diversification by vicariance as well as dispersal. Dispersal was also possible because the main part of the Galápagos lies in shallow water created by the broad Galápagos platform, which formed by the accumulation of lavas (^{Geist et al. 2008}).

Aegialomys galapagoensis, the Galápagos rice rat, is only found in one population on Santa Fé Island and is the latest species to have colonized the archipelago. Our results show high genetic diversity for the species and a signal of recent expansion. This finding is consistent with the fact that A. galapagoensis and N. narboroughi have been considered common by virtually all researchers visiting the islands (^{Clark 1980}; ^{Patton and Hafner 1983}; ^{Clark 1984}; ^{Key and Muñoz-Heredia 1994}). In fact, ^{Dowler et al. (2000)} reported that the species was abundant and rice rats could be seen running among the lava rocks before dark. Another reason is that both Santa Fé Island and Fernandina Island are the only two islands that do not have introduced species of rats and mice (^{Dowler et al. 2000}), allowing the persistence of the native species. This species is considered vulnerable by the Red List of the International Union for Conservation of Nature (^{IUCN 2019}). We dated the events of speciation and diversification to around 1.11 Ma and 111,200 years ago, respectively; both of them occurred after the emergence of Santa Fé Island, dated at 2.9 Ma (^{Geist et al. 2014}). ^{Clark (1980)} conducted an ecological study and found that A. galapagoensis has high survival and low reproduction relative to congeners of other geographic areas. He did not know if this strategy was specific to A. galapagoensis or a common feature of the Galápagos rodent life history (^{Harris and Macdonald 2007}).

Nesoryzomys narboroughi, the large Fernandina rice rat, is only found on Fernandina Island. It is the oldest species within the extant rodents of Galápagos, with a speciation event dated at 2.23 Ma; however, its diversification is the most recent among the extant Nesoryzomys, starting around 169,400 years ago. This species is considered vulnerable according to the Red List (^{IUCN 2019}); however, ^{Dowler et al. (2000)} verified the presence of sustainable populations, which are found from the coastline to the volcano rim. We found that the genetic diversity of this species is high, and their populations are under recent expansion. This could be due to the fact that Fernandina Island does not have introduced species of rats and mice, allowing the persistence of N. narboroughi.

Nesoryzomys fernandinae, the small Fernandina rice rat, is also found on Fernandina Island in sympatry with N. narboroughi (^{Hutterer and Hirsch 1980}; ^{Dowler et al. 2000}). Its populations range from the coast up to the volcano rim. ^{Dowler et al. (2000)} found that this species is more abundant at high elevations along the volcano rim, where vegetation, primarily Scalasia spp., is most dense. Meanwhile, N. narboroughi is more abundant near the coast, suggesting that the species distributions are influenced by habitat (^{Dowler and Carroll 1996}). N. fernandinae was considered extinct (^{IUCN 1996}; ^{Nowak 1999}), but now its status is vulnerable according to the Red List (^{IUCN 2019}). This species has the highest genetic diversity within the genus, and it shows a more stable population in comparison with N. narboroughi, which showed recent expansion. We dated the events of speciation and diversification to 1.58 Ma and 452,900 years ago, respectively. The presence of N. narboroughi and N. fernandinae on Fernandina Island suggest that volcanic activity has not occurred on an island-wide scale sufficient to destroy their populations as ^{Merlen (2014)} suggested. He also proposed that Fernandina Island is a refuge for biodiversity in Galápagos based on the formation of endemic species, the establishment of terrestrial ecosystems at several altitudes, and climate conditions.

Fernandina Island has the most active volcano and is the youngest major island in the archipelago. Its emergence is suggested to have occurred approximately 32,000 years ago with a maximum emergence around 60,000 years ago (^{Geist et al. 2014}). The dates that we obtained for the speciation and diversification of N. narboroughi and N. fernandinae are older than the dates proposed for the origin of the island.

Nesoryzomys swarthi, the Santiago Galápagos mouse, is endemic to Santiago Island where it exists as a single population in the arid zone of the north-central coast (^{Dowler et al. 2000}; ^{Harris et al. 2006}). This species coexists with Rattus rattus and Mus musculus on the island. However, the three species are only sympatric near the beach, where the densities of the invasive species are lower (^{Dowler et al. 2000}; ^{Harris et al. 2006}). ^{Brosset (1963)} considered that N. swarthi was extinct based on the introduction of diseases or parasites by invasive species, and/or competition with them. At present, it is considered vulnerable according to the Red List (^{IUCN 2019}). Our results showed that the species has relatively high genetic diversity, despite being the lowest genetic diversity among the other native rodents. Two of our analyses show a signal of recent expansion, but one of them rejected this hypothesis. We consider that more genetic analyses, including more loci, are necessary to resolve this discrepancy. It is possible that the species underwent a recent bottleneck; however, the genetic signal is not strong. A bottleneck event in the species would be consistent with the fact that it was considered extinct for decades while it went undetected until ^{Dowler et al. 2000} rediscovered it. However, we do not have enough data to thoroughly test this hypothesis. The emergence of Santiago Island has been dated to a minimum of 800,000 years ago and a maximum of 1.4 Ma, which agrees with the time of speciation and diversification that we found for N. swarthi around 1.58 Ma and 525,500 years ago, respectively. ^{Harris and Macdonald (2007)} performed an ecological study and found that the unpredictable environment of the Galápagos arid zone has selected for a strategy of high survivorship and low reproduction in N. swarthi that is more typically found among desert Heteromyidae than other members of the Oryzomyini. ^{Clark (1980)} found the same result in A. galapagoensis. It seems that this strategy is a common feature of native rodents of the Galápagos Islands.

According to our results, we have two main conclusions. The first is that populations of the extant rodents on the Galápagos show high genetic diversity and most show a signal of recent expansion. Despite the challenges that are faced by these species by the extreme arid conditions on the islands, the impact of human activities, and the presence of invasive species, our results suggest that these species do not have genetic signatures implying that they have experienced dramatic population declines and, rather, that populations are demographically stable. However, based on previous studies this conclusion should be considered carefully. For example, the desert-adapted life histories of these species are impressive because their ancestors were not well adapted to dry environments and oryzomyine rodents of semiarid habitats often lack the physiological adaptations to xeric conditions (^{Best 1988}; ^{Meserve 1978}; ^{Ribeiro et al. 2004}; ^{Harris and Macdonald 2007}). Moreover, the fact that the abundance of these endemic rodents is positively related to vegetation density suggests that the presence of these species is resource-limited (^{Clark 1980}; ^{Dowler and Carroll 1996}; ^{Harris and Macdonald 2007}). This could explain why some native species can coexist with invasive species or coexist in sympatry with other native species if there are enough resources to minimize competition. In this sense, there is now evidence that introduced rodents may not be the sole cause of decline in native rodents in Galápagos; feral cats may be a second important factor in rodent decline and extinction (^{Dowler et al. 2000}; ^{Dexter et al. 2004}), as well as loss of habitat and resources or the introduction of pathogens. In light of these previous findings and our results, we should continue to monitor these vulnerable species, performing more fine-scale ecological and population studies in the future.

Our second major conclusion is that the major speciation events of the four extant rice rats occurred within the archipelago during the Pleistocene. In the case of A. galapagoensis and N. swarthi, their diversification occurred on Santa Fé and Santiago Islands, respectively, and their speciation likely occurred within their current range. However, for the species inhabiting Fernandina Island, N. narboroughi and N. fernandinae, speciation and diversification likely occurred on a different island, because Fernandina Island is younger than these events. We hypothesize that after speciation on a different island, they dispersed to Fernandina Island. Isabela Island is the largest within the archipelago, and it is close to Fernandina. Its emergence is dated to around 500,000 to 800,000 years ago, and it is therefore possible that diversification of both species occurred there. However, we dated speciation events older than this emergence. Speciation of Nesoryzomys likely occurred on islands that emerged at least between 1.5 to 2.3 Ma. It is difficult to know where it occurred because it is recognized that volcanic islands leave little evidence of their ancient biological past, because the lava flows consume organic matter (^{Steadman 1986}). Thus, species may have gone through events of colonization and speciation and subsequently been driven to extinction, leaving no sign (^{Merlen 2014}). ^{Geist et al. (2014)} suggested that because the Galápagos archipelago is dynamic over evolutionary timescales, for any model of dispersal, colonization, speciation, and radiation involving island geography more than 20,000 years ago, the current map of the Galápagos Islands is close to irrelevant. Despite these problems, we suggest that future studies using genome-wide molecular markers and sampling of extinct species and fossils from the islands, as well as additional continental species, will add further insights into the origin and relationships of the endemic rodents of the Galápagos Islands, which remain ambiguous.