Pleistocene climatic oscillations in Neotropical open areas: Refuge isolation in the rodent Oxymycterus nasutus endemic to grasslands

Pleistocene climatic oscillations favoured the expansion of grassland ecosystems and open vegetation landscapes throughout the Neotropics, and influenced the evolutionary history of species adapted to such environments. In this study, we sampled populations of the rodent Oxymycterus nasutus endemic to open areas in the Pampas and Atlantic Forest biomes to assess the tempo and mode of population divergence using an integrative approach, including coalescence theory, ecological niche models, and morphometry. Our results indicated that these O. nasutus populations exhibited high levels of genetic structure. Six major mtDNA clades were found, structuring these biomes into distinct groups. Estimates of their divergence times was indicated to be 0.571 myr. The high degree of genetic structure is reflected in the analyses of geometric morphometric; skull differences between lineages in the two ecoregions were detected. During the last glacial maximum, there was a strong increase in suitable abiotic conditions for O. nasutus. Distinct molecular markers revealed a population expansion over time, with a possible demographic retraction during the post-glacial period. Considering that all clades coalesce with the last interglacial maximum, our results indicated that reduction in suitable conditions during this period may have resulted in a possible vicariance associated with refuge isolation.


Introduction
The glaciations occurred in the Pleistocene (2.58-0.01 million years ago [myr]), a period of climatic history that is known with near accuracy [1], which led to the migrations and reductions in the population sizes of several species, followed by recolonization and population expansion as glaciers retreated [2][3][4]. Therefore, the divergence promoted by glaciers have probably affected populations distributed in allopatry, which could be detected in genealogies, as observed in typical glacial refuges [5].

Sample collection
A total of 135 specimens of O. nasutus from 45 localities covering the entire known distribution of this species (Fig 1; Table 1) were surveyed. Field collection (permits granted by the Ministério do Meio Ambiente, Brazil-SISBIO 30204-1, 52650-1 and 29358-1) followed the guidelines of the American Society of Mammalogists for the use of wild mammals in research (https://www.mammalsociety.org/committees/animal-care-and-use). Specimens were euthanized via overdoses of isofluorane and cervical dislocation, procedures authorized by The Animal Care and Use Committee of the Institute of Biological Sciences at UFRGS, Brazil, for this study Specimens trapped were deposited both in the mammal collections of UFRGS and MCNU.

Phylogenetic relationships
DNA was isolated from 107 specimens using the PureLink Genomic DNA extraction kit (Invitrogen, Life Technologies), following the manufacturer's instructions. Polymerase chain reaction (PCR) was performed to amplify the mitochondrial DNA (mtDNA) cytochrome b (Cytb) gene (801 bp) and the nuclear locus beta-fibrinogen, intron 7 (Fgb-I7) (408 bp) according to the method described by Smith and Patton [24] and Matocq et al. [25]. The nuclear marker was chosen as it represents a single-copy locus, and is informative and evolving at a different rate compared to the mtDNA [26]. PCR products were stained with GelRed (Biotium) and checked on 1% agarose gel, purified with enzymatic method (Exonuclease and Alkaline Phosphatase; Amersham Biosciences, Piscataway, NJ) and sequenced with Sanger method using both primers (forward and reverse). All sequences obtained were deposited in GenBank under the accession numbers: Cyt b, MF766110 to MF766188 and Fgb-I7, MF766189 to MF766256.
Alignment and editing of the mtDNA Cytb sequences was performed in the Clustal W algorithm using the MEGA 7 software. For the nuclear Fgb-I7, sequences were aligned using MAFFT v.7.245 [26] with the auto setting. Variable sites of this nuclear marker were assessed in the original chromatograms to ensure correct identification of the heterozygotes. Heterozygous sites were identified when two different nucleotides were indicated at the same position in the chromatograms, with the weakest peak reaching at least 25% of the strongest signal. The IUPAC symbols were applied for coding the double peaks. The gametic phase of each haplotype was identified computationally using the program PHASE v2.1 [27,28] as implemented in DnaSP v5.10 [29] to reconstruct putative alleles of the nuclear marker for use in downstream haplotype-based analyses. The program was ran with 500 burn-in steps and 500 iterations, including the allowed intragenic recombination for the given data set. We also used a 0.6 output probability threshold for haplotypes and genotypes, as this was shown to reduce the number of unidentified genotypes with or without slight increase in the number of false positives [30].
The most appropriate substitution model for mtDNA phylogeny from the haplotypic data was the HKY+G that was determined based on the Akaike Information Criterion in the jMo-delTest v2.1.10 software [31]. The phylogeny was reconstructed using Bayesian Inference (BI) in the software BEAST 1.8.4 [32]. Two independent Markov Chain Monte Carlo (MCMC) runs, each with four streams per 25 million steps of the MCMC, sampled every 1000 generations, and discarding 2.5 million burn-in (about 10% trees discarded), starting the initial trees with randomness, without restriction were performed. Speciation tree prior was modelled on Birth-Death Process using the BEAST software. Parameter convergence was checked in Tracer v1.5 [33] as to whether effective sample sizes (ESS) reached 200. The remaining trees were used to calculate the posterior probabilities for each node.
The burn-in was determined in Tracer v1.6 based on the trajectory parameters, and 10% of the initial trees were removed and summarized in TreeAnnotator. The consensus tree generated was visualized and edited in Figtree v1.4.3 (http://tree.bio.ed.ac.uk/software/figtree/). There were no fossils for Oxymycterus that could be used for a specific estimate of the substitution rate. Hence, we used an evolutionary rate already published [34] of 2.37% per million years (Myr -1 ) for the genus. Accordingly, the divergence times were estimated in BEAST assuming an uncorrelated relaxed clock model and a normal distribution for the substitution rate, with a mean of 2.37% Myr-1, and a standard deviation of 0.25% Myr-1, to allow some uncertainty in the evolutionary rate. One representative sequence for each haplotype was used. Sequences of Oxymycterus delator (U03525), Oxymycterus dasytrichus (AF454768), Oxymycterus amazonicus (AF454765), and sister lineages of O. nasutus [35], were included in the analysis. The time of the most recent common ancestor (TMRCA) for relevant nodes and major mitochondrial clades was reported as the mean value of node height with 95% highest posterior density interval (HPD). The nodes were supported with the posterior probability (PP) of ! 90% [36]. We also reconstructed the phylogenetic tree among nuclear DNA (nDNA) haplotypes (Fgb-I7) data set with phased alleles and single-copy sequences using BI in BEAST. However, due to the low variability and lack of support, we performed downstream analysis based on the haplogroups inferred by mtDNA. The topology obtained from Fgb-I7 evolutionary tree is shown in S1 Fig. The evolutionary relationships between the haplotypes for each data set (i.e. mtDNA and nDNA) were estimated using the median-joining method implemented in the Network v5.0 [37] (http://www.fluxus-engineering.com). For assessing the evolutionary relationships of the Fgb-I7 gene, haplotypes were identified by coalescent-based Bayesian method implemented in PHASE v2.1. Missing/gaps sites were verified and invariant sites were removed from the dataset. Indels were treated as single mutational events.

Intraspecific variation and historical demography
Statistical parameters of diversity, such as the number of variable sites (S), number of haplotypes (h), haplotype diversity (Hd), nucleotide diversity (π), and average number of nucleotide differences (k) were measured using DNASP v5.10 [29]). We calculated the genetic distance between the recognized haplogroups of Cytb (average distances) using p-distance model with 1000 bootstrap replications using MEGA v7 [38]. The defined clades were assumed based on the phylogeny obtained through Bayesian analysis.
Historical demographic changes such as signatures of demographic expansion, and equilibrium or decline were examined for the species and for each haplogroup using neutrality tests, such as Tajima's D [39] and Fu's Fs [40] statistics, employing the software DNASP v5.10. Pairwise mismatch distribution was performed in DNASP v.5.10 to infer the historical demography of O. nasutus, calculated with the expected frequency based on a population growthdecline model. The sum of squared deviations (SSDs) between the observed and expected mismatch distribution and the raggedness index (r) were calculated to test the null hypothesis of spatial expansion using ARLEQUIN v3.5 [41]. In addition, we performed a Bayesian skyline plot (BSP) analysis, which does not assume a priori any growth model and infers the effective population size through time based on coalescent theory [42]. The BSP was used to estimate the dynamics of alterations in the overall population size over time for the O. nasutus dataset. BEAST v1.8.4 [32] was used for estimating the BSP of mtDNA as described above, with minor alterations in the number of parameters sampled at every 10,000 steps. For nDNA, a posterior distribution model of effective population size through time was generated using a MCMC sampling scheme. Two independent analyses were ran for 40 x 10 7 generations (sampling at every 10,000 and 10% burn in) under a HKY+G substitution model (obtained through the jModelTest 2 software), assuming a relaxed molecular clock model. However, we used the substitution rates and dates estimated from previous studies to calibrate a relaxed molecular clock and approximate divergence times for the main phylogroups retrieved in our study. We implemented a lognormal prior distribution based on the substitution rate for the cricetid genus Peromyscus (0.006 ± 0.003 substitution/site/million years) for the Fgb-I7 data set [43]. Skyline reconstruction was performed in Tracer v1.5, and the median and 95% credibility interval were plotted as a function of time.

Ancestral area reconstruction
Aiming to reconstruct the ancestral range of O. nasutus, we performed a Statistical dispersalvicariance analysis (S-DIVA) on the maximum clade credibility tree generated from the BEAST analysis using the software RASP 3.02 [44][45][46]. The DIVA method [44] develops an approach based on parsimony events and reconstructs ancestral distributions based on a simple biogeographic model and a three-dimensional cost matrix. Penalties are not assigned when speciation is the result of vicariance; however, dispersal and extinction events have a penalty of one per unit area added to or deleted from a distribution. S-DIVA method rectified the problems in DIVA analysis and suggested possible ancestral ranges at each node and also calculated the probabilities for each ancestral range at the nodes. Due to the current area of occurrence and based on historical aspects of the region (expansion and retreat of forest formations), we reconstructed the ancestral areas only within the Pampas (A) and Atlantic Forest/grassland mosaic (B) regions. According to the results obtained by BI and the presence of restricted haplogroups in a single ecoregion, the clades recognized in O. nasutus were categorized in these areas. We used the output files generated by BEAST (total and condensed trees) for the analysis.

Morphological analyses
Digital photographs were taken from the dorsal, lateral, and ventral views of the skull of 89 specimens of O. nasutus (Table 1). Only adult specimens were photographed to minimize the ontogenetic effects (the complete eruption of the third molar was the criterion to separate juveniles from adults). Two-dimensional digital images were taken using a Canon PowerShot G10 camera with 14.7 megapixels resolution (4416 x 3312) in the macro function of the automatic mode and without flash or zoom. The pictures were taken from a standard distance of 127 mm for all the specimens. A total of 17 landmarks were digitized in the dorsal view, 16 landmarks in the lateral view, and 32 landmarks in the ventral view (S2 Fig), using the TpsDig2 software [47]. The same person (WTP) conducted the digitization. The choice of landmarks was based on previous studies with sigmodontines [48,49]. In the dorsal and ventral view, all individuals were marked on both sides of skull. The description of all the landmarks employed is given in S1 Appendix. After digitization, a Generalized Procrustes Analysis (GPA) was conducted on the matrix of landmark coordinates to remove the effects of scale, position, and orientation [50]. The centroid size, derived by the square root of the sum of squared distances of each landmark from the centroid of the configuration [51], was used as a measure of size. We appended the natural log-transformed centroid size into the matrix of shape coordinates to work in the form of space for further analyses.
Our goal was to investigate whether morphology is divergent among the haplogroups found in the genetic analyses (Northwest, Central, Eastern, Steppes Plain, Southern, and Taim Wetland), and also between physiognomies or "environmental groups" (Pampas versus Atlantic Forest). First, we conducted a principal component analysis (PCA) in the form matrix for each view. The number of PCs necessary to achieve 100% variation in each view (16 PCs for dorsal view, 31 PCs for ventral view, and 29 PCs for lateral view) were used as response variables for downstream analyses. We performed multivariate analyses of variance (MANOVA) and discriminant analyses with Jackknife cross validation (DA) to investigate how morphology was structured among the six genetic haplogroups (first predictor) and the two environmental groups (second predictor). MANOVAs and DAs were performed independently for each predictor and for each skull view. All procedures were carried out in the software R [52]), with the packages geomorph [53] and Morpho [54]. Visualization of shape changes was made by comparing the shapes among groups using discriminant functions implemented in MorphoJ [55]. Geographic distribution maps and modelling Ecological niche modelling (ENM) was carried out in MAXENT v 3.3.3e [56] to predict suitable present and past potential distribution areas of O. nasutus based on climatic variables. Such modelling has been performed favourably compared to other analytical alternatives for presence-only data [56][57][58][59]. Information on the geographic distribution of the species was based on 51 occurrence records obtained from literature and collections (S2 Appendix). We verified the accuracy of the geo-referenced data in ArcGIS v 10.0. To avoid data clustering [60], we limited our database to a single record per km 2 . For modelling settings, the function 'Auto features' was selected, and the distributions were modelled through the 'cross-validate' parameter, applying a maximum number of iterations at 500. We verified model performance using the area under the 'Receiver Operating Characteristic (ROC) Curve' (AUC) calculated by MAXENT. Values between 0.7 and 0.9 indicated good discrimination [61].
Projections for the past climatic conditions were developed for three periods: LGM c. 21 ka; Last IinterGlacial (LIG) c. 120-140 ka; and mid-Holocene c. 6 ka. Current and past distributions were modelled using 19 bioclimatic data layers available from the WorldClim database (http://www.worldclim.org) at 30 arc-sec resolution (~1 km 2 ), with the exception of the layers

Intraspecific genetic variation and evolutionary patterns
We identified 73 variable sites in the mtDNA dataset resulting in 38 haplotypes and 32 variable sites in the Fgb-I7 locus resulting in 40 haplotypes (with one 2-bp indel) (Tables 1 and 2). Standard diversity indices (haplotype diversity, nucleotide diversity, mean number of pairwise differences) for both the markers are presented in Table 2.
Bayesian consensus tree based on mtDNA dataset depicted six highly supported (BPP > 0.90) haplogroups corresponding to two distinct biomes: Northwest, Central, and Eastern clades are distributed along the Atlantic Forest biome, and Steppes Plain, Taim Wetland, and Southern clades assigned to the Pampas biome. However, the relationship between these groups was not entirely clear, considering that they did not cluster in major clades exclusive to each of the two biomes (Fig 2). The haplotype network based on the mtDNA also supported such tree structure. Nonetheless, despite the several mutational steps between haplogroups, the haplotype network showed several median vectors indicating non-sampled or extinct ancestral sequences (Fig 3). The two main widespread haplogroups were Central and Southern clades. The Central clade occurred throughout the domain of Atlantic Forest, mainly dispersed over the Araucaria forests and the mosaic of forest/highland grasslands (Fig 2). This clade was  (Table 1). Only Piraquara locality in Parana state, Brazil, comprised haplotypes from two distinct mtDNA clades (Fig 1). nDNA network revealed 36 low-frequency haplotypic variants, with a reticulate evolutionary relationship. The majority of haplotypes differed by one substitution site. None of the six major clades in the mitochondrial tree was recovered for the Fgb-I7 sequences. However, nuclear haplotype H1 was present in almost all mtDNA clades (except Taim Wetland). All mtDNA clades showed the presence of unique alleles, but Taim Wetland group shared the haplotype 31 with Steppes Plain and Southern clades.
The genetic distances between the Bayesian clades ranged from 1% to 2.5% for mtDNA, whereas it was close to 0 for Fgb-I7 (S1 Table).
The most common ancestor for all clades of O. nasutus was estimated to the middle Pleistocene (0.5715 myr; 95% HPDs = 0.3657-0.8471 mya) (Fig 2, Table 3). The mtDNA haplotypes (n = 38) clustered into six lineages; individual clades diverged between 265.8 and 147.3 myr (HPDs = 0.438.8-0.046.7 myr). The ancestral haplotype could not be inferred for this marker due to the presence of several median vectors in the network. However, the BEAST-derived tree indicated H7 (Northwest clade) as the oldest clade. Tajima's D and Fu's Fs tests were negative and non-significant for both mtDNA and nDNA markers, indicating that O. nasutus might have experienced a recent population expansion ( Table 2).    Table 2). The results of the mismatch distribution for mtDNA and nDNA analysis was approximately unimodal (Fig 4). Non-significant SSD statistic (SSD = 0.00112542, P = 0.944) and raggedness index value (r: 0.0032, P = 0.9770) under the spatial expansion models failed to reject the spatial expansion model. Similar patterns suggested demographic expansion for nDNA (SSD = 0.00322295, P = 0.570; r: 0.01631205, P = 0.780). Strong evidence of demographic expansion came from the BSP. The results showed different patterns between the mitochondrial and the nuclear datasets. For the mtDNA, a constant population size was noted over the last 0.375-1 mya after a long phase of demographic stability; the population then appeared to have experienced an accelerated demographic expansion phase approximately 0.030-0.300 myr, followed by a decrease in population size after 0.030 myr. For the nDNA marker, results showed a slow growth in the effective population size across time (0.400-1.2 mya), demographic stability for 30-400 myr, with smooth decrease posteriorly (Fig 4).  skull (S3 Fig). The percentage of correct classification among the haplogroups were: 62.92% for the dorsal view, 65.16% for the ventral view, and 55.05% for the lateral view. Among the environmental groups, the percentages of correct classification were: 93.25% for the dorsal view, 82.02% for the ventral view, and 80.89% for the lateral view. The percentages of correct classification were similar between the genetic and the environmental haplogroups, considering the total number of groups in each one (i.e. six haplogroups would result in~16% of the correct classification by chance, while with two groups this percentage will rise to 50%; deviations from the random classification are similar between the predictors). Detailed results of discriminant analyses can be found in the S2 and S3 Tables.  ENM results (Fig 6) and the potential LGM distributions were more widespread. On the other hand, the potential distribution of MID-HOL was reduced, which showed more similarities with the current conditions. The potential distributions of the LIG were the most restricted of all the ENMs. Suitable habitats in the LGM showed a great expansion in relation to the LIG. However, in all the models tested, areas of low suitability present between the Pampas and

Phylogeography of O. nasutus
Atlantic Forest domains were identified, as well as in the western part of the territory of Uruguay. In addition, the coastal areas did not have suitable conditions during the LGM.
According to S-DIVA analysis (Fig 7), Pampas ecoregion (A/AB, node 1) was the most likely region for the origin of the dispersal of O. nasutus. Dispersal events probably progressed from Pampas domain (A) to areas that are currently translocated into the Atlantic Forest domains (B) and also ranging to other regions of the Pampas. Vicariance signals, detected on node 2 (AB), indicated that vicariance events occurred between the Northwest clade and the other haplogroups that had dispersed between the Pampas and the Atlantic Forest. A new dispersion event was detected on node 3 (A/AB), scattering the haplogroups (Southern, Steppes Plain, Central, and Eastern) along the Pampas and the current Atlantic Forest areas. Finally, on node 4, a new vicariance was verified between the haplogroups that were present on the Atlantic Forest areas (Central and Eastern) and the Pampas (Steppes Plain).

Genetic and morphological variation in the Pampas and Atlantic Forest
The absence of Cytb haplotypes that was common to Atlantic Forest and Pampas groups was intriguing. However, the scenario was intricate, as a lack of further structuring between the two biomes in the Bayesian and median-joining network analysis was evident. This pattern indicated a possible isolation following the historical scenario of gene flow, as depicted by the Fgb-I7 network. Patterns of skull form were also roughly in agreement with the molecular findings. The differentiation between the physiognomies of O. nasutus populations were evident in the classification after the cross-validation. Subtle differences in the skull shape of O. nasutus populations were also depicted among the environments. Such skull differences might be related to genetic grouping or with environmental factors, exerting some form of local selection and/or guiding phenotypic plasticity in each physiognomy [8]. Similarly, studies on haplogroups of sympatric akodonts (e.g., Deltamys kempi [63,64] and Scapteromys spp. [65,66]). were also divergent in the geometric or linear morphometric analysis of the skull.
A conspicuous geomorphological feature of the connection among the southern Atlantic Forest and Pampas landscapes is the escarpment of Serra Geral (Meridional Plateau), which might have acted as a long-term barrier to gene flow between the Pampas and the current Atlantic Forest open areas at the eastern point of contact of these two biomes. Thus, the altitudinal gradient formed might have influenced certain level of differentiation, thereby rendering a group of lineages dispersing across the highlands of the Meridional Plateau, and others through the sedimentary lowlands of the Paraná basin, and later spreading along the coastal plain. Mountain ranges and ridges have been considered as topographic elements associated with speciation processes in akodontine rodents [67][68][69][70][71]. Since the differentiation in O. nasutus populations is recent (time-tree revealed that the diversification probably began in the Middle Pleistocene) and the divergence process seems incomplete, it can be considered at the population-level and not at species-level yet.
The population structure of O. nasutus was quite complex; therefore, it is possible that the recent expansion of forest formations of the South Atlantic Forest, such as the Araucaria Forest and the seasonally dry Forest lato sensu [21], along with the escarpment of Serra Geral have created a new geographical barrier that separated the Atlantic Forest and Pampas lineages. Hence, divergence in the mtDNA clades suggested that the lineages have persisted in isolation across time, characterizing an 'grassland refuge' pattern. Accordingly, the species adapted to dry and/or cold environments and restricted to small areas surrounded by unsuitable habitats can persist through interglacial microrefugium [72].

Phylogeographic patterns
The phylogenetic tree revealed a more complex evolutionary history in the open areas of Pampas and Atlantic Forest instead of simple reciprocally monophyletic groups. Six major mtDNA clades were observed, which fall within the range of intraspecific divergence described for the Akodontini tribe [63,64,71,73]. The mtDNA haplotype network analysis in this study was congruent with such tree structure, representing the six haplogroups separated by several mutational steps, and inserted in specific ecoregions. However, median joining vectors present in the network did not rule out the hypothesis of non-sampled or possibly extinct haplotypes. Although differentiation was less evident, patterns found in nDNA haplotypes provided indication of a demographic expansion event. Similarly, results of mtDNA analysis also suggested demographic growth. In addition, intraspecific Fgb-I7 variability was conspicuously lower than Cytb variability, as reported in other studies [43,63,74]. Molecular markers such as mtDNA and nDNA generally show different phylogeographic patterns for the same biogeographical history owing to differentiation in effective population sizes, recombination, and mutation rates [75,76].
The estimates of genetic distances for the mtDNA clades were moderately high (up to 2%). This might be the result of an ongoing process of isolation, and because mtDNA accumulates substitutions 5-10 times faster than the single-copy nuclear DNA [75]. The highest distances were observed among Taim Wetland versus other clades and between Southern and Northwest clades; however, none of these genetic distances were recovered for the nuclear locus (divergence close to zero). The highest genetic distance observed for these clades (Taim Wetland, Northwest, and South) might relate to the ancestral condition of these groups, as depicted by BI and mtDNA network, along with the fact that they are geographically distant from each other. According to the S-DIVA analysis, tandem dispersion events followed by vicariance might have probably influenced the phylogenetic pattern of mtDNA observed in the study.
The Northwest clade was the northernmost lineage; so, the effects of the dynamics of forest formations/open areas might have reached earlier in this group, making it more suitable to environmental pressure. Positive values for the Tajima's D test (0.44358, P = 0.7550) for this clade (from mtDNA) and other negative but non-significant values (for both markers), indicated a possible contraction effect on this lineage. The divergence between the Steppes Plain and Central and Eastern clades dates back to 0.350 myr (0.220-0.505 myr). Moreover, some forest formation similar to Seasonal Forest, which at present is adjacent to the escarpment of the Serra Geral might have acted concomitantly. This hypothesis was supported by the results of the ENM that indicated an area of very low suitability between the zone of contact of these two ecoregions. This is consistent with the escarpment of Serra Geral, which might have corroborated mainly due to the vicariance between the Steppes Plain and Central and Eastern clades, as detected by S-DIVA analyses. Similarly, a strong phylogenetic break in this region was evident for other akodonts recently described, such as Scapteromys meridionalis [71] and Deltamys araucaria [64]. An evident pattern is the genetic divergence within the Southern clade in the Pampas ecoregion (Uruguay). In this study, ENM supported up to a limited dispersion, indicating a possible historical barrier acting in this region. In this case, the Rio Negro could be possibly acting as a barrier for Southern clade, wherein areas of low suitability were identified in all the prediction models. Southern clade was one of the most spread throughout the Pampas ecoregion. Although neutrality test indicated expansion, Southern clade was the only clade where suitable areas were identified mainly in the coastal regions and continental shelf, according to the LGM niche models.
A remarkable phylogeographic aspect was the genetic divergence of the Taim Wetland clade. The mtDNA dataset revealed a statistically supported phylogenetic break between the Taim Wetland clade and other clades, indicating a possible historical barrier acting in this region. This lineage originated about 0.147 myr (95% HPD = 0.047-0.290 myr) and was the first to have dispersed in the S-DIVA analysis. This clade was inserted on the plains of southernmost Brazil (RS), where climatic oscillations occurred during the middle Pleistocene and Holocene, resulting in marine transgressions and regressions of relative sea level that shaped the South Atlantic Coastal Plain (SACP). The formation of the SACP is related to the sedimentary processes associated with the sea transgressive events known as 'Barriers', which occurred in the middle Pleistocene with the formation of Barrier I (~0.400 myr), Barrier II (~0.325 myr), Barrier III (~0.120 myr), and Holocene, when Barrier IV (~0.005 myr) was established [77]. Possibly, this clade accessed this region after the formation of barriers II and III, between 0.325 and 0.120 myr, respectively. The geological evolution of the "Barriers" also formed the Patos-Mirim complex, the largest lagoon-barrier system in South America [77]. In addition, three paleochannels related to the Jacuí, Camaquã, and Jaguarão rivers were identified in the middle and southern coastal zone of RS [78]. These paleochannels were associated with the large persistent hydrographic elements such as Mirim Lagoon, São Gonçalo Channel, and the Patos Lagoon estuarine channel, which might have possibly contributed to the isolation of this lineage in the southernmost domains of the SACP. Similarly, recent studies show a strong phylogenetic break in the same region for Scapteromys tumidus [65] and Deltamys kempi [64]. Thus, major elements of Patos-Mirim lagunar complex might constitute a geographical barrier for the historical gene flow.

Climate change and landscape
Overall, substantial demographic changes associated with Pleistocene climatic oscillations were observed in more than half of the biota investigated in South America. Considering only the taxa to be associated with open vegetation, a total of 68% of the studied species experienced population expansion during the glacial periods [6]. The historical biogeographic approach presented herein appeared to corroborate this pattern detected for the open areas-dwelling species.
All divergence times within O. nasutus populations estimated by BEAST occurred during the Middle Pleistocene, wherein a minimum of eight glaciations occurred in the Middle-Late Pliocene in the southern Andes of South America [79]. However, the greatest glaciations occurred during the early Pleistocene, such as the Great Patagonian Glaciations, between 1.16 and 1.01 mya, when the glaciers advanced up to 200 km east of the Andes mountains, and stretching along the Pacific and south of Atlantic coast [79,80]. After the Great Patagonian Glaciations, 13 minor glacial and interglacial periods were recorded for this region in the Early-Middle Pleistocene that reached a maximum around 25,000 and 16,000 years BP (see [79]).
Although the emergence of O. nasutus occurred in the earlier periods (ca. 500 kyr), the patterns showed by the present and past niche models depicted expansion of high probability areas during the LGM (~21,000 years BP), favouring its distribution. Probably, a pattern similar to LGM might have occurred during the Great Patagonian Glaciations, where the effects of glaciations might have increased the open areas and retracted the forest formations. The areas of Southern Atlantic Forest were dominated by grasslands between 42,000-10,000 years BP; the period comprising the LGM. Forest elements were restricted to sites in deep river valleys and in coastal lowlands, indicating a cold and dry climate [21]. In general, grasslands adapted to cold conditions prevailed in the southern and southeastern Brazil until around 11,500 years BP [81]. Since O. nasutus is intrinsically related to open areas and is abundantly found in these habitats [17][18][19], cyclic events of the dynamic expansion and retraction of open areas and forests that occurred during glacial and interglacial periods have probably caused species dispersal and historical distribution. Such expansions of open areas promoted by glaciations made the environment more suitable for establishment of new areas, which was detected by S-DIVA analysis. It was interesting to reveal that no suitable areas in the coastal zone were increased during LGM according to the niche model. Palynological studies on the coasts of the states of Santa Catarina and Paraná indicated that fields were abundant in this region and also on the exposed continental shelf during the LGM period, whereas tropical tree species were practically absent [23]. The continental shelf during LGM presented a dry climate, hampering the colonization by any species adapted to humid/wet habitats. Areas of moderate suitability were found in the coastal region of southern Brazil and Uruguay during the LIG, Middle Holocene, and current models, with a reduction of suitability in the inland areas. Thus, the coastal regions in the Pampas ecoregion might have served as a refuge for O. nasutus during the interglacial periods.
O. nasutus is known to inhabit the subtropical grassland in wet or flooded seasons [17]. Accordingly, the species is adapted to humid environments, being found in bunch grass and in tall grass near the streams and rivers [82], sandbanks near wet lands across coastal semi-fossorial habit [19], and coastal vegetation in the South Atlantic Coastal Plain [83]. Thus, in addition to grassland, O. nasutus might have followed river routes in search of humid environments. Such habitat specialization might reflect the dispersion ability to more extensive regions during the glacial periods.
In this context, we suggest that the two major rivers in Rio Grande do Sul state (the Jacuí and Ibicuí) might have acted (and still be acting) as river barrier for the Steppes Plain clade. The haplotypes of this clade are distributed across the western edge of the Patos Lagoon and the interior of the plains of southern Brazil. Therefore, ENM suggested that these areas might be unsuitable for the restriction caused mainly by the river Jacuí. Others rivers could have affected the limits of dispersion in O. nasutus. Rivers Uruguay and Paraná also might have acted in limiting the dispersion during the glacial periods. The paleoecological niche modelling suggested that at LGM, the River Uruguay limited the dispersion for Steppes Plain and Southern clades in Uruguay and Southern Brazil to reach new areas to the west. Similarly, the River Paraná limited the dispersion for Northwest and Central Clade in highlands. Finally, Paraná, and Paranapanema rivers might have acted as physical barriers in the dispersion to the north, mainly for the Northwest clade, due to their amplitude and water volume. These inferences were based on the ENM results, wherein areas of low suitability were verified in all the prediction models.