https://orcid.org/0000-0002-4473-7640
Abstract
The Pampas of South America represents one of the most extended subtropical–temperate grasslands in the world. This ecoregion was influenced by Pleistocene climate oscillations. The sigmodontine rodent assemblage inhabiting this region is a good model system to analyse the impact of these climate changes on evolutionary histories. We performed a comparative phylogeographical study on seven species to evaluate the concordance of phylogeographical patterns among them, to assess if they experienced congruent and synchronous demographic changes, and posit putative centres of origin and dispersal routes. Four species (Calomys musculinus, Oligoryzomys flavescens, Oxymycterus nasutus and Oligoryzomys nigripes) showed evidence of demographic expansion. In the first three, population increases started during Marine Isotope Stage 5 (MIS 5) while in Ol. nigripes, the rise started during MIS 7; this rise would have continued to the present in all cases. Necromys lasiurus, Scapteromys tumidus and Scapteromys aquaticus did not show a pattern of expansion. Four centres of origin were identified; in general, populations sharing the same putative refugium followed common dispersal routes. Our results indicate that the Pampas offered relatively stable habitats over time, mainly in areas associated with watercourses or coastlines, suggesting that recent Pleistocene climate oscillations have had a moderate impact on this ecoregion compared to other regions of South America.
INTRODUCTION
The geographical distribution of a taxon is the result of the interaction among geological, environmental and ecological processes in space and time. Therefore, studies of the geographical distribution of faunal units, i.e. an assemblage of local fauna that shares ecological attributes and has experienced the same events over a period of time (Petronio & Marcolini, 2013), should integrate all those processes. This approach can be addressed by comparative phylogeography, which studies demographic processes and evolutionary divergence patterns of multiple taxa across a common region to evaluate if the species have responded to shared historical events in a similar way (Edwards et al., 2022). Theoretically, co-distributed species exposed to common environmental and geological changes would have congruent phylogeographical patterns. However, those patterns can also be affected by species’ ecological characteristics such as dispersal abilities and food or habitat preferences, and/or by the possible instability over time of the assemblage itself, which would conceal the effects of their common history (Taberlet et al., 1998; Zhang et al., 2012).
In South America, the landscape configuration has been greatly influenced by Pleistocene climate oscillations between cold and generally dry (glacial) and warm and generally humid (interglacial) periods, combined with marine transgressions that shifted the shorelines, giving rise to complex scenarios of species diversification (Colinvaux et al., 1996; Behling, 2002; Antonelli & Sanmartin, 2011). Despite this, Beheregaray (2008) highlighted the scarcity of phylogeographical research in South America. Most of such studies refer to species from tropical regions and, to a lesser extent, to those from temperate regions. In the subtropical regions of South America, the role of geological events and climate oscillations on species demographic processes and evolutionary divergence patterns have not been analysed in depth (Turchetto-Zolet et al., 2013). One of such regions is the Pampas, one of the most extended subtropical–temperate grasslands in the world, encompassing more than 750 000 km2 in Uruguay, central–eastern Argentina and southern Brazil (Bilenca & Miñarro, 2004; see Fig. 1). The dominant physiognomy of this ecoregion consists of herbaceous plants, although other vegetation types such as gallery forests, shrublands, wetlands or meadows are also represented. Based on this diversity, several authors have distinguished different subregions within the Pampas (Soriano et al., 1992; Olson et al., 2001). However, in the present study, we consider the Pampas as a whole as in Cabrera & Willink (1973). Regarding the effects of Pleistocene climatic oscillations, it has been suggested that this ecoregion stretched over forested areas during glacial periods and was invaded by forest biomes during interglacial periods, thus reducing its geographical range (Behling, 2002; Lorenz-Lemke et al., 2010; Tonello & Prieto, 2010). In addition, the lowering of sea level during glacial periods exposed additional landmasses on the Atlantic coast, also favouring Pampas expansion (Clapperton, 1993; Tomazelli & Villwock, 1996, 2000).
Sampling sites of rodent specimens (coloured dots) in the Pampas ecoregion (grey).
Phylogeographical studies in this ecoregion have focused mainly on single species (e.g. the plant Turnera sidoides, Moreno et al., 2018; the ant Acromyrmex striatus, Cristiano et al., 2016; the frog Pseudopaludicola falcipes, Langone et al., 2016), with very few covering multiple taxa (e.g. plants of the genera Petunia and Calibrachoa, Fregonezi et al., 2013). The assemblage of sigmodontine rodents that inhabit the Pampas is one of the best known in South America in relation to its composition, competitive interactions, and species ecological characteristics such as habitat preferences, diet and daily activities (Busch & Kravetz, 1992; Bilenca & Kravetz, 1995; Pardiñas et al., 2004). This assemblage mainly comprises species of the genera Akodon, Calomys, Necromys, Oligoryzomys, Oxymycterus, Reithrodon and Scapteromys, which are primarily considered omnivorous, consuming overall plants and invertebrates, with high trophic overlap among them (Bilenca et al., 1992; Ellis et al., 1998). However, they also show seasonal variations in their diets which have been attributed to mechanisms of coexistence in the assemblage (Castellarini et al., 1998). Furthermore, by analysing the distribution of congeneric species, Bilenca (1993) found a predominance of only one of these species at the same time and in the same location, which was also considered a mechanism to reduce interspecific competition, allowing their co-occurrence in the Pampas ecoregion. Based on this knowledge and in the rodents’ shared morphophysiological adaptations (i.e. size, rostrum elongation and/or mandible robustness and masticatory apparatus), Pardiñas et al. (2010) proposed that the structure of the pampean assemblage has been conditioned by a common evolutionary history. On the other hand, the species within this assemblage differ in habitat preferences, some being endemic to open areas, others being more associated with open gallery forests adapted to flooded areas, and others being extreme generalists in their habitat type (Patton et al., 2015). Their ecological similarities, together with, on the one hand, the suggested coexistence mechanisms among them, and, on the other, their ecological differences, make the rodent assemblage from the South American Pampas a good model to analyse the impact of Pleistocene climatic oscillations in the region through a comparative phylogeographical approach. Therefore, we aimed to (1) analyse the phylogeographical patterns in rodent species inhabiting the Pampas, (2) assess if they experienced congruent demographic changes, (3) estimate the time during which these changes occurred and whether they were synchronic among the species and (4) posit putative shared centres of origin and dispersal routes among these rodents. Our prediction is that the pampean rodent assemblage will exhibit concordant phylogeographical patterns and demographic history, as a result of the influence of shared climatic and ecological factors on the evolution of their populations.
MATERIAL AND METHODS
Sample collection
We used the mitochondrial cytochrome b gene (Cyt-b) as a molecular marker due to its high mutation rate at the intraspecific level (Bradley & Baker, 2001). It is also the most widely used gene in phylogeographical studies of rodents, which makes it suitable for comparative analyses. We used data from previous research from our laboratory in Calomys musculinus (González-Ittig et al., 2007), Oligoryzomys flavescens (Rivera et al., 2018) and Oligoryzomys nigripes (González-Ittig et al., 2010; Colombo et al., 2018). In addition, new sequences were obtained from tissues of specimens of C. musculinus and Necromys lasiurus following DNA extraction, amplification and sequencing as in González-Ittig et al. (2010). Other available Cyt-b sequences of these and other species belonging to the pampean sigmodontine assemblage were obtained from the GenBank database. Therefore, Oxymycterus nasutus, Scapteromys aquaticus and Scapteromys tumidus were also included since they had adequate sample size data to perform phylogeographical analyses. Other rodent species inhabiting the Pampas could not be analysed due to their low sample size (Akodon azarae, Calomys laucha and Oxymycterus rufus) or because we were not able to correlate the haplotypes and localities from the original article (Quintela et al., 2017) with the sequences published in GenBank (Deltamys kempi). Subterranean rodents of the genus Ctenomys were also excluded, since their distribution and demographic processes are influenced mainly by edaphic conditions (Kittlein & Gaggiotti, 2008; Mora et al., 2013) rather than by climatic oscillations during the Pleistocene. Several of the studied species have a distribution that extends beyond the Pampas ecoregion; however, in accordance with our aims, only sequences from rodents captured in Pampas localities or in others nearby were considered. Details of species, capture sites, geographical coordinates, museums or collections where voucher specimens are deposited, specimen identification codes or voucher numbers and GenBank accession numbers are given in Supporting Information Table S1.
Studied rodents
Calomys musculinus (common name: drylands vesper mouse) has a wide geographical distribution in Argentina, Bolivia and western Paraguay (Supporting Information, Fig. S1). This rodent prefers microhabitats with low herbaceous vegetation and croplands (Salazar-Bravo, 2015), and is one of the most abundant rodents in the pampean agroecosystems. It is considered an opportunistic species, showing a great capacity for adaptation to anthropogenic environments (Busch et al., 2000).
Necromys lasiurus (common name: hairy-tailed Bolo mouse) is distributed in lowlands of tropical, subtropical and temperate regions of central and eastern South America (Pardiñas et al., 2015; Supporting Information, Fig. S1). It is a cursorial, mainly diurnal and habitat generalist species (Vieira et al., 2005; Becker et al., 2007).
Oligoryzomys flavescens (common name: yellow pygmy rice rat) is a species complex distributed in Argentina, Bolivia, Uruguay, Paraguay, and south and eastern Brazil (Rivera et al., 2018; Supporting Information, Fig. S1). In this study, we analyse only specimens that correspond to the Ol. flavescens flavescens cluster, which is distributed mostly in the Pampas and Espinal (Rivera et al., 2018). This rodent has terrestrial habits, preferring areas close to streams and humid lowlands with dense vegetation, and stable environments with tall weeds in disturbed habitats (Mills et al., 1991; Busch & Kravetz, 1992; Bonaventura et al., 2003).
Oligoryzomys nigripes (common name: black-footed pygmy rice rat) is found from central and eastern Brazil and eastern Paraguay to Uruguay and north-eastern and eastern Argentina (Weksler & Bonvicino, 2005; Cirignoli et al., 2006; Supporting Information, Fig. S1). It prefers areas with abundance of shrubs, linked to its scansorial behaviour (Vieira & Monteiro-Filho, 2003).
Oxymycterus nasutus (common name: long-nosed hocicudo) is endemic to grasslands in Uruguay and south-eastern Brazil (Supporting Information, Fig. S1). Its most common habitats are wet grasslands, coastal sand banks, steppes and plateaus (Oliveira & Goncalves, 2015). It exhibits a diurnal pattern of activity (Paise & Vieira, 2006).
Scapteromys aquaticus (common name: Argentine swamp rat) is found in central–eastern and north-eastern Argentina, southern Paraguay, and on the coast of the Uruguay river in western Uruguay and south-western Brazil (D’Elia & Pardiñas, 2004; Bonvicino et al., 2013; Supporting Information, Fig. S1). Scapteromys tumidus (common name: swamp rat) ranges across most of Uruguay and along the eastern coastal strip in Brazil (Freitas et al., 1984; D’Elía & Pardiñas, 2004; Fig. S1). Both species prefer semi-aquatic habitats and areas near watercourses associated with flooded microhabitats with low plant coverage (Bonaventura et al., 2003). These rodents display a nocturnal pattern of activity (Hershkovitz, 1966; D’Elía & Pardiñas 2004).
Phylogeographical analyses
For each species, multiple sequence alignments were performed with Muscle (Edgar, 2004) using default parameters. Mitochondrial haplotypes were defined with DNAsp v.6.12.03 (Rozas et al., 2017); positions from the DNA sequences containing missing data were excluded from the analysis. The number of sequences (n), number of haplotypes (NHap), matrix length (ML), average number of nucleotide differences (K), nucleotide diversity (π) and haplotype diversity (Hd) were calculated for each species with Arlequin v.3.5.1.3 (Excoffier & Lischer, 2010). Haplotype networks were constructed using the Median-Joining algorithm implemented in PopART v.1.7 (Leigh & Bryant, 2015). To assess the past demographic history of the species, Tajima’s D and Fu’s FS neutrality tests and mismatch distributions (Harpending, 1994) were also calculated with Arlequin v.3.5.1.3. Bayesian skyline plot (BSP) analysis, computed with Beast v.1.10 (Suchard et al., 2018), was used to infer past population demography. A sequence evolution model for each species was selected using jModeltest v.2.1.3 (Darriba et al., 2012) following the Bayesian information criterion (BIC). An uncorrelated relaxed clock was used with an initial substitution rate of 2.3%/ Myr estimated by Smith & Patton (1993) for the Cyt-b gene for South American rodents. Calibrations were implemented in the form of lognormal prior distributions as suggested by Ho & Phillips (2009) and the following log(Mean) (μ), log(Stdev) (δ) and offset (this is a hard minimum bound, or when available, the age of the fossils): Ol. flavescens μ = −0.8, δ = 0.6, offset = 0.126; Ol. nigripes μ = −0.5, δ = 0.8, offset = 0.126; and N. lasiurus (the fossil was originally under Necromys benefactus name) μ = −0.5, δ = 0.6, offset = 0.126. With these values, the 95% tail of the prior distribution falls around 2 Myr to create a soft maximum bound (for details of the fossils see Pardiñas et al., 2002). For the rest of the species, fossil data were not available. The model parameters were included in Beast, and the analysis was run for 100 million generations. Chain convergence, effective population size estimations, parameter confidence intervals and BSP reconstructions were checked using Tracer v.1.7 (Rambaut et al., 2018). A Mantel correlation test between pairwise FST and geographical distances among localities was performed to search for a pattern of isolation by distance. This test was performed, after excluding localities with only one sequence, with 10 000 random permutations in zt v.1.1 (Bonnet & Van de Peer, 2002).
To reconstruct the spatio-temporal history of lineage dispersal, we used a relaxed random walk model (RRW), implemented in Beast v.1.10 (Lemey et al., 2010). The priors for sequence evolution were the same as described for the BSP analysis. We used the Coalescent GMRF Bayesian Skyride model as a tree prior (Minin et al., 2008). Runs were performed with 600–800 million generations depending on the convergence of the parameters in the data set of each species. The convergence of a parameter was evaluated by its effective sample size (ESS) calculated in Tracer v.1.7. A maximum clade credibility tree was generated with the software TreeAnnotator of the Beast package, with a burnin of 10% of the trees. Spatio-temporal reconstructions were obtained using Spread v.1.0.7 (Bielejec et al., 2016). The kml file output with the phylogeographical reconstruction was visualized in Google Earth Pro (Google Inc.).
RESULTS
Sample collection
A total of 30 sequences of C. musculinus, 27 of N. lasiurus, 63 of Ol. flavescens, 24 of Ol. nigripes, 76 of Ox. nasutus, 16 of S. aquaticus and 128 of S. tumidus captured in localities of the Pampas or in other nearby localities were analysed (Fig. 1, Supporting Information, Table S1; Fig. S1). Given their differences in length, sequences were trimmed to build matrices in a range between 600 and 800 bp. Matrix length, number of haplotypes, and nucleotide and haplotype diversity for each species are given in Table 1.
Sample information and genetic diversity indices for the Cyt-b gene of seven rodent species of the pampean assemblage.
| Species | N | NHap | ML (bp) | K ± SD | S | π ± SD | Hd ± SD |
|---|---|---|---|---|---|---|---|
| Calomys musculinus | 30 | 17 | 600 | 2.5 ± 1.4 | 13 | 0.004 ± 0.003 | 0.94 ± 0.02 |
| Necromys lasiurus | 27 | 9 | 801 | 3.6 ± 1.9 | 15 | 0.004 ± 0.004 | 0.68 ± 0.09 |
| Oligoryzomys flavescens | 63 | 30 | 711 | 3.8 ± 1.9 | 35 | 0.005 ± 0.003 | 0.93 ± 0.02 |
| Oligoryzomys nigripes | 24 | 17 | 801 | 7.3 ± 3.5 | 35 | 0.009 ± 0.001 | 0.97 ± 0.02 |
| Oxymycterus nasutus | 76 | 31 | 738 | 7.9 ± 3.2 | 49 | 0.009 ± 0.006 | 0.94 ± 0.02 |
| Scapteromys aquaticus | 16 | 5 | 801 | 2.9 ± 1.6 | 9 | 0.004 ± 0.002 | 0.77 ± 0.07 |
| Scapteromys tumidus | 128 | 11 | 789 | 4.3 ± 2.2 | 16 | 0.006 ± 0.003 | 0.85 ± 0.01 |
| Species | N | NHap | ML (bp) | K ± SD | S | π ± SD | Hd ± SD |
|---|---|---|---|---|---|---|---|
| Calomys musculinus | 30 | 17 | 600 | 2.5 ± 1.4 | 13 | 0.004 ± 0.003 | 0.94 ± 0.02 |
| Necromys lasiurus | 27 | 9 | 801 | 3.6 ± 1.9 | 15 | 0.004 ± 0.004 | 0.68 ± 0.09 |
| Oligoryzomys flavescens | 63 | 30 | 711 | 3.8 ± 1.9 | 35 | 0.005 ± 0.003 | 0.93 ± 0.02 |
| Oligoryzomys nigripes | 24 | 17 | 801 | 7.3 ± 3.5 | 35 | 0.009 ± 0.001 | 0.97 ± 0.02 |
| Oxymycterus nasutus | 76 | 31 | 738 | 7.9 ± 3.2 | 49 | 0.009 ± 0.006 | 0.94 ± 0.02 |
| Scapteromys aquaticus | 16 | 5 | 801 | 2.9 ± 1.6 | 9 | 0.004 ± 0.002 | 0.77 ± 0.07 |
| Scapteromys tumidus | 128 | 11 | 789 | 4.3 ± 2.2 | 16 | 0.006 ± 0.003 | 0.85 ± 0.01 |
N: number of sequences; NHap: number of haplotypes; ML: matrix length; K: average number of nucleotide differences; S: number of variable sites; π: nucleotide diversity; Hd: haplotype diversity; SD: standard deviation.
Sample information and genetic diversity indices for the Cyt-b gene of seven rodent species of the pampean assemblage.
| Species | N | NHap | ML (bp) | K ± SD | S | π ± SD | Hd ± SD |
|---|---|---|---|---|---|---|---|
| Calomys musculinus | 30 | 17 | 600 | 2.5 ± 1.4 | 13 | 0.004 ± 0.003 | 0.94 ± 0.02 |
| Necromys lasiurus | 27 | 9 | 801 | 3.6 ± 1.9 | 15 | 0.004 ± 0.004 | 0.68 ± 0.09 |
| Oligoryzomys flavescens | 63 | 30 | 711 | 3.8 ± 1.9 | 35 | 0.005 ± 0.003 | 0.93 ± 0.02 |
| Oligoryzomys nigripes | 24 | 17 | 801 | 7.3 ± 3.5 | 35 | 0.009 ± 0.001 | 0.97 ± 0.02 |
| Oxymycterus nasutus | 76 | 31 | 738 | 7.9 ± 3.2 | 49 | 0.009 ± 0.006 | 0.94 ± 0.02 |
| Scapteromys aquaticus | 16 | 5 | 801 | 2.9 ± 1.6 | 9 | 0.004 ± 0.002 | 0.77 ± 0.07 |
| Scapteromys tumidus | 128 | 11 | 789 | 4.3 ± 2.2 | 16 | 0.006 ± 0.003 | 0.85 ± 0.01 |
| Species | N | NHap | ML (bp) | K ± SD | S | π ± SD | Hd ± SD |
|---|---|---|---|---|---|---|---|
| Calomys musculinus | 30 | 17 | 600 | 2.5 ± 1.4 | 13 | 0.004 ± 0.003 | 0.94 ± 0.02 |
| Necromys lasiurus | 27 | 9 | 801 | 3.6 ± 1.9 | 15 | 0.004 ± 0.004 | 0.68 ± 0.09 |
| Oligoryzomys flavescens | 63 | 30 | 711 | 3.8 ± 1.9 | 35 | 0.005 ± 0.003 | 0.93 ± 0.02 |
| Oligoryzomys nigripes | 24 | 17 | 801 | 7.3 ± 3.5 | 35 | 0.009 ± 0.001 | 0.97 ± 0.02 |
| Oxymycterus nasutus | 76 | 31 | 738 | 7.9 ± 3.2 | 49 | 0.009 ± 0.006 | 0.94 ± 0.02 |
| Scapteromys aquaticus | 16 | 5 | 801 | 2.9 ± 1.6 | 9 | 0.004 ± 0.002 | 0.77 ± 0.07 |
| Scapteromys tumidus | 128 | 11 | 789 | 4.3 ± 2.2 | 16 | 0.006 ± 0.003 | 0.85 ± 0.01 |
N: number of sequences; NHap: number of haplotypes; ML: matrix length; K: average number of nucleotide differences; S: number of variable sites; π: nucleotide diversity; Hd: haplotype diversity; SD: standard deviation.
Phylogeographical analyses
Haplotype networks are shown in Figure 2. Calomys musculinus exhibited three haplotypes in central positions connected through a few mutations (in most cases only one) to several exclusive haplotypes (Fig. 2A). Oxymycterus nasutus showed a complex haplotype network presenting three haplogroups separated by several mutation steps, one distributed in the centre and north of the sampling area (haplogroup North), another mostly in the south (haplogroup South) and the last one between the other two (Fig. 2B). Haplogroup North showed a star-like pattern with one frequent, widely distributed haplotype and several others connected, in general, by a few mutational steps; haplogroup South did not show any defined geographical structure. The third haplogroup did not include enough haplotypes to infer a network structure. Oligoryzomys flavescens showed a two-star-like pattern with the majority of the haplotypes connected by one to three mutational steps; one of the central haplotypes had a widespread distribution and the other one was exclusive to one population (Fig. 2C). The Ol. nigripes network showed two geographical haplogroups separated by several mutational steps: one that clustered haplotypes from Argentina (localities 1 and 2, haplogroup West) and the other, those from Brazil and Uruguay (haplogroup East) (Fig. 2D). Within haplogroup East there was a complex pattern in which haplotypes were separated by several mutations; however, a star-like pattern could be inferred around the haplotype present in localities 5 and 10 (Fig. 2D). Scapteromys tumidus, S. aquaticus and N. lasiurus did not show a star-like pattern (Fig. 2E, F and G, respectively). The S. tumidus network exhibited several frequent haplotypes, all of them having wide geographical distributions. In the S. aquaticus network, most of the haplotypes were exclusive to one locality, diverging among them by a few mutational steps. Necromys lasiurus had one very frequent haplotype widely distributed in a non-central position, separated by one to four mutations from the other haplotypes.
Median-joining networks showing the relationships and relative abundance of cytochrome b haplotypes of: A, Calomys musculinus; B, Oxymycterus nasutus; C, Oligoryzomys flavescens; D, Oligoryzomys nigripes; E, Scapteromys tumidus; F, Scapteromys aquaticus; and G, Necromys lasiurus. Each bar through the solid line represents a 1 nucleotide difference between haplotypes. Empty rhombuses represent median vectors, that is hypothetical intermediate haplotypes (unsampled or extinct). For each species, numbers correspond to the sampling localities indicated in Supporting Information Table S1. Haplogroups in Oxymycterus nasutus and Oligoryzomys nigripes are framed with a dotted grey line.
Demographic analyses
Pronounced genetic and geographical divergence among haplogroups may indicate different evolutionary histories. Hence, the following analyses were performed separately for haplogroup East in Ol. nigripes and haplogroups North and South in Ox. nasutus, which included enough sequences to infer demographic processes (Fig. 2E and B, respectively). The remaining sequences of these species were not included in the analyses. For the rest of the species, all sequences were considered. In C. musculinus, Ol. flavescens and Ol. nigripes haplogroup East, and in haplogroups North and South of Ox. nasutus, Fu’s FS neutrality index was negative and significant, suggesting population expansion. Tajima’s D neutrality index was negative and significant only for Ol. flavescens and Ol. nigripes haplogroup East. For the rest, Tajima’s D index was also negative but non-significant (Table 2). In N. lasiurus, S. aquaticus and S. tumidus, neutrality tests were in general positive but non-significant, indicating that populations of these species are in equilibrium. Mismatch distribution analyses showed a classical unimodal curve only for C. musculinus. The rest of the species and haplogroups had multimodal or irregular curves (Supporting Information, Fig. S2), although Rg (Harpending’s raggedness index) and SSD (deviations from the sum of squares) were all non-significant (Table 2). Hence, the null hypothesis of population expansion was not rejected.
Neutrality tests and statistics of the mismatch distribution analysis for the Cyt-b gene of the seven rodent species representing the pampean assemblage.
| Species | D | FS | SSD | Rg | MSE |
|---|---|---|---|---|---|
| Calomys musculinus | −0.77 (NS) | −11.51* | 0.010 (NS) | 0.056 (NS) | HKY+I |
| Necromys lasiurus | −0.07 (NS) | 1.29 (NS) | 0.109 (NS) | 0.191 (NS) | HKY |
| Oligoryzomys flavescens | −1.60* | −19.96* | 0.005 (NS) | 0.014 (NS) | HKY+I |
| Oligoryzomys nigripes haplogroup East | −1.58* | −5.10* | 0.004 (NS) | 0.014 (NS) | HKY |
| Oxymycterus nasutus haplogroup North | −1.37 (NS) | −5.77* | 0.003 (NS) | 0.013 (NS) | HKY+I |
| Oxymycterus nasutus haplogroup South | 1.02 (NS) | −8.025* | 0.003 (NS) | 0.015 (NS) | HKY |
| Scapteromys aquaticus | 0.19 (NS) | 1.28 (NS) | 0.043 (NS) | 0.087 (NS) | HKY |
| Scapteromys tumidus | 1.28 (NS) | 2.71 (NS) | 0.010 (NS) | 0.022 (NS) | HKY |
| Species | D | FS | SSD | Rg | MSE |
|---|---|---|---|---|---|
| Calomys musculinus | −0.77 (NS) | −11.51* | 0.010 (NS) | 0.056 (NS) | HKY+I |
| Necromys lasiurus | −0.07 (NS) | 1.29 (NS) | 0.109 (NS) | 0.191 (NS) | HKY |
| Oligoryzomys flavescens | −1.60* | −19.96* | 0.005 (NS) | 0.014 (NS) | HKY+I |
| Oligoryzomys nigripes haplogroup East | −1.58* | −5.10* | 0.004 (NS) | 0.014 (NS) | HKY |
| Oxymycterus nasutus haplogroup North | −1.37 (NS) | −5.77* | 0.003 (NS) | 0.013 (NS) | HKY+I |
| Oxymycterus nasutus haplogroup South | 1.02 (NS) | −8.025* | 0.003 (NS) | 0.015 (NS) | HKY |
| Scapteromys aquaticus | 0.19 (NS) | 1.28 (NS) | 0.043 (NS) | 0.087 (NS) | HKY |
| Scapteromys tumidus | 1.28 (NS) | 2.71 (NS) | 0.010 (NS) | 0.022 (NS) | HKY |
D: Tajima’s neutrality index; FS: Fu’s neutrality index; SSD: deviations from the sum of squares; Rg: Harpending’s raggedness index statistic; MSE: model of sequence evolution. Asterisks indicate significant values at P < 0.05; NS: non-significant.
Neutrality tests and statistics of the mismatch distribution analysis for the Cyt-b gene of the seven rodent species representing the pampean assemblage.
| Species | D | FS | SSD | Rg | MSE |
|---|---|---|---|---|---|
| Calomys musculinus | −0.77 (NS) | −11.51* | 0.010 (NS) | 0.056 (NS) | HKY+I |
| Necromys lasiurus | −0.07 (NS) | 1.29 (NS) | 0.109 (NS) | 0.191 (NS) | HKY |
| Oligoryzomys flavescens | −1.60* | −19.96* | 0.005 (NS) | 0.014 (NS) | HKY+I |
| Oligoryzomys nigripes haplogroup East | −1.58* | −5.10* | 0.004 (NS) | 0.014 (NS) | HKY |
| Oxymycterus nasutus haplogroup North | −1.37 (NS) | −5.77* | 0.003 (NS) | 0.013 (NS) | HKY+I |
| Oxymycterus nasutus haplogroup South | 1.02 (NS) | −8.025* | 0.003 (NS) | 0.015 (NS) | HKY |
| Scapteromys aquaticus | 0.19 (NS) | 1.28 (NS) | 0.043 (NS) | 0.087 (NS) | HKY |
| Scapteromys tumidus | 1.28 (NS) | 2.71 (NS) | 0.010 (NS) | 0.022 (NS) | HKY |
| Species | D | FS | SSD | Rg | MSE |
|---|---|---|---|---|---|
| Calomys musculinus | −0.77 (NS) | −11.51* | 0.010 (NS) | 0.056 (NS) | HKY+I |
| Necromys lasiurus | −0.07 (NS) | 1.29 (NS) | 0.109 (NS) | 0.191 (NS) | HKY |
| Oligoryzomys flavescens | −1.60* | −19.96* | 0.005 (NS) | 0.014 (NS) | HKY+I |
| Oligoryzomys nigripes haplogroup East | −1.58* | −5.10* | 0.004 (NS) | 0.014 (NS) | HKY |
| Oxymycterus nasutus haplogroup North | −1.37 (NS) | −5.77* | 0.003 (NS) | 0.013 (NS) | HKY+I |
| Oxymycterus nasutus haplogroup South | 1.02 (NS) | −8.025* | 0.003 (NS) | 0.015 (NS) | HKY |
| Scapteromys aquaticus | 0.19 (NS) | 1.28 (NS) | 0.043 (NS) | 0.087 (NS) | HKY |
| Scapteromys tumidus | 1.28 (NS) | 2.71 (NS) | 0.010 (NS) | 0.022 (NS) | HKY |
D: Tajima’s neutrality index; FS: Fu’s neutrality index; SSD: deviations from the sum of squares; Rg: Harpending’s raggedness index statistic; MSE: model of sequence evolution. Asterisks indicate significant values at P < 0.05; NS: non-significant.
BSP analyses showed an increase in the effective population size in four of the seven species analysed (C. musculinus, Ol. flavescens, Ol. nigripes haplogroup East and Ox. nasutus haplogroup North). For N. lasiurus, S. aquaticus and S. tumidus and also for Ox. nasutus haplogroup South there was no evidence of demographic change (Supporting Information, Fig. S3). The best-fitting model of sequence evolution for this analysis is shown in Table 2. To evaluate the temporal congruence of the demographic changes, we plot together the four lineages that experienced an effective population size increase (Fig. 3). To establish a correlation between demographic changes and Pleistocene climate oscillations, we added the time series of globally distributed benthic δ18 O records (Lisiecki & Raymo 2005) into Figure 3. Two different patterns could be seen: a population increase that started about 120 000 years ago during Marine Isotope Stage 5 (MIS 5) and has continued until the present in C. musculinus, Ol. flavescens and Ox. nasutus haplogroup North, while in Ol. nigripes haplogroup East, the increment started about 250 000 years ago during MIS 7. Both MIS 5 and MIS 7 correspond to interglacial periods.
Combined Bayesian Skyline Plots showing mean estimates (thick lines) and 95% confidence intervals (dashed lines) of effective population size through time for the four rodent species that showed demographic expansion. The x-axis represents time in thousand years before the present (Kya). The y-axis represents the log conversion of the effective population size. Green lines: Calomys musculinus; orange lines: Oligoryzomys flavescens; blue lines: Oligoryzomys nigripes haplogroup East; and red lines: Oxymycterus nasutus haplogroup North. Top panel: inverted time series of globally distributed benthic δ18 O records (Lisiecki & Raymo 2005) with marine isotope stages (MIS) indicated. Blue shading indicates glacial periods.
Mantel tests could not be performed in N. lasiurus, Ox. nasutus haplogroup South or S. aquaticus, since few localities remained when those with one individual were excluded from pairwise FST comparisons. The five species analysed showed a significant correlation between FST and geographical distances (C. musculinus, r = 0.85; P < 0.05; Ol. flavescens, r = 0.45; P < 0.05; Ol. nigripes haplogroup East, r = 0.73; P < 0.05; Ox. nasutus haplogroup North, r = 0.56; P < 0.05; S. tumidus, r = 0.53; P < 0.05) (Supporting Information, Fig. S4).
Lineage dispersal
The phylogeographical reconstruction of lineage dispersal showed that Ol. nigripes haplogroup East originated 460 kya (thousand years ago) on the north-west coast of Patos lagoon and later dispersed to the south of the lagoon 250 kya, reaching its southern end about 130 kya. Inland spread to the north-west would have occurred 65 kya and only recently to the south-west (Fig. 4A). Oxymycterus nasutus haplogroup North originated around 530 kya also on the west coast of Patos lagoon and then expanded inland to the south about 440 kya. It would then have extended northwards and southwards from 400 to 30 kya, following the west coast of the lagoon. Finally, dispersal from the southern inland populations to the west and from the west coast of Patos lagoon populations to the north occurred ~140 and 60 kya, respectively (Fig. 4B). On the other hand, Ox. nasutus haplogroup South would have originated in the south-eastern sea coast of Uruguay 390 kya and later migrated slightly to the south 250 kya and extended inland to the north 160 kya. The southern populations would have dispersed to the north-east from 150 kya to recent times, thus colonizing the eastern sea coast of Uruguay (Fig. 4A). Scapteromys tumidus populations would also have arisen on the south-eastern sea coast of Uruguay, although more recently than the rest of the studied rodent species (nearly 90 kya). They would have spread to the north-east 80 kya (to the south of Patos lagoon) and to the south and south-east, following the coastline, from 80 to 50 kya. Then, the southern populations might have migrated to the north-east, east and north 60, 45 and 40 kya, respectively. Finally, multiple short dispersal events would have occurred among the coastal populations 30 kya, leading to the genetic structure observed today (Fig. 4B). Scapteromys aquaticus populations appeared in the lower Paraná river area ~720 kya, although the phylogeographical reconstruction of lineage dispersal shows migration events only recently, towards the north, south and west (Fig. 4B). The lower Paraná river area would also have been the centre of origin of the Ol. flavescens populations, but 300 kya. These populations may have dispersed westwards close to 210 kya; later, they reached the mouth of de la Plata river in the south nearly 170 kya and migrated to the north about 100 kya. The western populations might then have experienced a new dispersion further to the west, from 40 kya to the present. The lower Paraná river populations would have continued expanding to the east, south and south-west up to recent times (Fig. 4A). Calomys musculinus populations would have originated in central–eastern Argentina 245 kya and then colonized the south-eastern Pampas 160 kya, the Paraná river area in the north 140 kya and an inland area towards the north-west 120 kya. The southern populations might have migrated to the east following the coast 110 kya and the northern populations to the west about 75 kya (Fig. 4B). Necromys lasiurus populations would also have originated in central–eastern Argentina 580 kya and migrated to the south sea coast of Buenos Aires province 470 kya and to the north, near the Paraná river coast, 390 kya. Finally, the southern populations expanded to the west, reaching the westernmost point of the sampled distribution ~40 kya (Fig. 4A).
Location of the centres of origin (refugia) and lineage spatial diffusion retrieved with the Bayesian spatio-temporal diffusion analysis for each rodent species. Ellipses show the estimated centres of origin for each species or clade within the species, which are represented in different colours (one for each lineage). Arrows indicate the direction of lineage diffusion, and the estimated timing of arrival to the area is indicated in thousand years before the present. Arrows without numbers represent current dispersion events. dlPR: de la Plata river, ML: Mirim lagoon, NR: Negro river, PL: Patos lagoon, PR: Paraná river, SR: Salado river, UR: Uruguay river.
DISCUSSION
Demographic processes and the influence of Pleistocene climatic oscillations
To evaluate the impact of Pleistocene climate oscillations on the evolutionary history of the assemblage of South American Pampas rodents, we performed a comparative phylogeographical analysis among seven species from the region. Previous phylogeographical studies had been carried out in some of these species but encompassing a smaller (C. musculinus, González-Ittig et al. 2007), the same (Ol. flavescens, Rivera et al. 2018; S. tumidus, Quintela et al. 2015), or a broader (Ol. nigripes, Miranda, 2007; Ox. nasutus, Peçanha et al., 2017) geographical area than the one considered here. No phylogeographical studies had previously been performed for N. lasiurus or S. aquaticus from the region.
The results showed that C. musculinus, Ol. flavescens, Ol. nigripes haplogroup East and Ox. nasutus haplogroup North presented mitochondrial star-like haplotype networks (Fig. 2), indicating a relatively recent demographic expansion in the Pampas followed by population differentiation by drift, with currently low to moderate levels of gene flow. These four lineages also showed a negative and significant Fu’s FS index, although Tajima’s D index was significant only for two of them (Ol. flavescens and Ol. nigripes haplogroup East). Calomys musculinus was the only species that presented a typical unimodal mismatch distribution; however, those of Ol. flavescens, Ol. nigripes haplogroup East and Ox. nasutus haplogroup North would also be compatible with an expansion model (Supporting Information, Fig. S2). BSP analyses also reveal an increase in the effective population size for these taxa (Fig. 3).
In contrast, the networks of N. lasiurus, Ox. nasutus haplogroup South, S. aquaticus and S. tumidus are compatible with a non-expansion pattern in which historically low levels of gene flow predominated. The neutrality tests were non-significant, except for Fu’s FS in Ox. nasutus haplogroup South. The curve obtained in the mismatch distribution analysis for this lineage, although irregular, could also be interpreted as indicating an pattern of expansion (Supporting Information, Fig. S2). By contrast, the multimodal curves obtained for N. lasiurus, S. aquaticus and S. tumidus are in line with a stability model. BSP analyses indicated that these rodents have maintained their effective population size throughout the Middle Pleistocene to the present (Fig. S3).
In general, our results are in agreement with those of the previous phylogeographical studies quoted above, with the exception of for S. tumidus. Based on a BSP analysis, Quintela et al. (2015) suggested that this species experienced a population expansion, although the neutrality tests were non-significant. Differences in parameter settings could have led to the discrepancy between BSP results. Hence, we performed the BSP analysis with parameter settings following Quintela et al. (2015), but again no expansion signal was observed (results not shown). In our study, positive and non-significant neutrality tests (Table 1), a clearly multimodal curve obtained in the mismatch distribution analysis (Supporting Information, Fig. S2), and a constant population size over time in the BSP analysis (Fig. S3) did not indicate a population expansion for S. tumidus in the Pampas ecoregion.
An pattern of isolation by distance was observed in the five species in which the Mantel tests could be performed. Scatterplots corresponding to C. musculinus, Ol. flavescens, Ol. nigripes haplogroup East and Ox. nasutus haplogroup North showed the typical pattern of regional equilibrium between gene flow and drift (Supporting Information, Fig. S4). However, in S. tumidus two groups of points are clearly shown: at small geographical distances, low and non-significant values of genetic differentiation were observed for some population pairs, while at all geographical distances, maximum and significant values of genetic differentiation were found; this pattern is in line with the dominance of drift in the species. In S. aquaticus the Mantel test was not performed because of the small sample size. However, almost all populations presented a different dominant and exclusive haplotype, even at small geographical distances, indicating restricted gene flow among them. Our results in the two Scapteromys species could be attributed to the ecological preference of swamp rats for flooded microhabitats. In a 2-year population ecology study of S. aquaticus (originally named S. tumidus), Bonaventura et al. (2003) found important changes in its abundance associated with a dry spring–summer period. In hydrophilic species such as these, the impact of precipitation regimes on population parameters could be strong, generating successive local extinctions followed by short-distance colonizations from few individuals. RRW analyses showed time estimations of the dispersal routes of S. aquaticus and S. tumidus in the order of tens of thousand years, while the rest of the species presented time estimations in the order of hundreds of thousand years (Fig. 4), which reinforces the idea that ecological characteristics of these species would explain their phylogeographical patterns rather than the environmental and geological history of the region.
Regarding the time of expansion, BSP analyses indicated a pronounced demographic increase phase beginning about 120 kya (125 000–100 000) in C. musculinus, Ol. flavescens and Ox. nasutus haplogroup North during the last maximum interglacial at MIS 5, and which continued throughout the glacial MIS 4, interstadial MIS 3, the Last Glacial Maximum (LGM, MIS 2) and the Holocene, up to the present. In Ol. nigripes haplogroup East, the analysis showed a slow growth in the effective population size over time, from 250 kya at MIS 7 throughout all the glacial and interglacial stages to the present (Fig. 3). This result suggests that although these four species underwent the same demographic process, the population expansions were asynchronous, starting at different interglacial periods. The reconstruction of population size over time was previously estimated in Ol. flavescens (Rivera et al. 2018) with similar results, and in S. tumidus (Quintela et al., 2015) for which we did not detect any expansion. In Ox. nasutus, Peçanha et al. (2017) found six haplogroups highly divergent in populations from the Pampas and Atlantic Rain Forest, which were analysed altogether, and a demographic expansion at about 300 kya was estimated. Here, we studied two of these six haplogroups, finding a population increase at about 120 kya for Ox. nasutus haplogroup North (Steppes plain lineage in Peçanha et al., 2017) and demographic stability in Ox. nasutus haplogroup South (Southern lineage in Peçanha et al., 2017). If the analysed population comprises a pool of different lineages, the demographic inference is likely to show a mixed pattern that reflects their diverse evolutionary histories (Malhi et al., 2001), which would explain the difference between the two studies. The evolutionary independence of the demographic patterns among clades or subspecies with different geographical distributions in the Pampas and, therefore, exposed to distinct historical environmental conditions, has been observed in other phylogeographical studies in plants (Turchetto et al., 2014; Moreno et al., 2018), a frog (Barrasso, 2014) and a gecko (Felappi et al., 2015), among others.
Several studies have associated species demographic changes in the Pampas with Pleistocene climatic oscillations (for a review see Turchetto-Zolet et al., 2013), although few of them have performed estimations of population size over time. Those few studies showed highly discordant results regarding when the demographic changes among species took place. Two different trends can be recognized regarding the effects of Pleistocene climatic oscillations in the Pampas. On the one hand, demographic increases by the end of the Pleistocene or during the Holocene have been reported in several species. For example, subterranean rodents of the genus Ctenomys showed demographic stability until the end of the LGM (~24 kya) in C. torquatus (Roratto et al., 2014) or until the Early Holocene (~8 kya) in C. talarum (Mora et al., 2013), when important increases in the effective population size occurred in the south of the Pampas. Several plant species from the Pampas in Uruguay and southern Brazil presented a similar historical pattern. Petunia integrifolia subsp. depauperata showed a smooth but constant population growth since ~20 kya in the LGM (Ramos-Fregonezi et al., 2015); Calibrachoa heterophylla and Petunia axillaris axillaris would have expanded during the Pleistocene–Holocene transition at ~12 kya (Mäder et al., 2013) and 10 kya (Turchetto et al., 2014), respectively. In general, the authors suggested that these expansions could be explained by the occurrence of novel environments free of competitors, such as sand dune or coastal habitats resulting from cycles of marine transgression/regression. Later, this process would also have been favoured by the increase in temperature and humidity during the Holocene.
On the other hand, the effective population size of several species grew or remained stable from the Middle Pleistocene to the present. Among species showing the latter is the neotropical angiosperm Eugenia uniflora, which experienced small or no changes in effective population sizes throughout the Pleistocene in riparian forests of the Pampas ecoregion (Turchetto-Zolet et al., 2016). The authors suggested that E. uniflora has been stable throughout the last 1 Myr in this region due to the survival of fragmented populations in multiple refugia even during glaciation peaks when the climate was dry and cold. In our study, S. aquaticus, S. tumidus, N. lasiurus and Ox. nasutus haplogroup South showed historical stability. However, due to their distinct ecological characteristics and non-overlapping distributions, it is not possible to infer a common scenario to explain their phylogeographical pattern.
Among the species that started their expansion in the Pleistocene around MIS 7, as estimated here for Ol. nigripes haplogroup East, we note the crab Aegla uruguayana that inhabits the bottom of streams, rivers and freshwater lakes in the Pampas. In this species, the clade that corresponds to A. uruguayana sensu stricto experienced demographic growth over the last 300 kyr (Zimmermann et al., 2021). In the gecko Homonota uruguayensis, the central clade within its distribution range (at the boundary region between Uruguay and Brazil) expanded ~250 kya (Felappi et al., 2015). For the frog Physalaemus henselii, Barrasso (2014) estimated a demographic expansion 250 kya for the clade distributed in Uruguay. Langone et al. (2016) observed that the P. falcipes clade from Brazil remained stable until about 250 kya, after which it experienced exponential growth. They suggested that Pleistocene glacial epochs had little impact on P. falcipes populations, as the species continued expanding with time. The same explanation can be extended to Ol. nigripes haplogroup East. Due to its scansorial behaviour, Ol. nigripes is closely associated with riparian forests of the Pampas. Similarly, the species A. uruguayana, P. henselii and P. falcipes are distributed predominantly in humid habitats. Hence, it is possible that climatic oscillations from the Middle Pleistocene to the present did not significantly affect environments near water bodies.
Several Pampas species have expanded more recently but prior to the LGM, such as C. musculinus, Ol. flavescens and Ox. nasutus haplogroup North. In the grasshopper Dichroplus elongatus from the eastern margin of the Paraná river a remarkable increase in population size occurred at ~150 kya (Rosetti & Remis, 2012). The plant Petunia axilaris parodii experienced a population growth at 100 kya (Turchetto et al., 2014) as did the perennial herb T. sidoides that showed a gradual increase in effective population size from about 75 kya after a long period of demographic stability (Moreno et al., 2018).
In most Pampas species, including those analysed here, no severe population bottlenecks were reported; on the contrary, there is evidence of population stability throughout the Pleistocene or of demographic expansions spanning glacial and interglacial cycles. This reinforces the idea that climatic oscillations, at least from the Middle Pleistocene to the present, did not have a strong impact on the species of this region. The lower latitude of the Pampas compared with temperate regions such as Patagonia could explain a diminished effect of glacial periods in the former. In addition, the Pampas is less continental and hence there is an increased maritime influence leading to lower temperature fluctuations and less drought in comparison with tropical Amazonia (Carnaval et al., 2009; Fitzpatrick et al., 2009). A mild impact of Pleistocene climatic oscillations in the Pampas would also explain the discordant and/or asynchronous demographic patterns, not only among rodents in our present study but also among other species in the ecoregion.
Centres of origin and dispersal routes
The RRW indicated four common centres of origin for the populations of the species analysed here: the west coastal plains of Patos lagoon in Brazil, the south-east coast of Uruguay, the lower Paraná river and the Salado river basin in Argentina, each one shared by two species. These areas would have acted as refugia for the species during unfavourable climate periods, from where populations later dispersed following specific routes.
The centre of origin in the west coastal plains of Patos lagoon was shared by Ol. nigripes haplogroup East and Ox. nasutus haplogroup North, both dated at ~500 kya (Fig. 4). This region has been reported as the centre of origin for populations of the frogs P. henselii (Barrasso, 2014) and P. falcipes (Langone et al. (2016), and of the plant T. sidoides (Moreno et al., 2018). Furthermore, Silva-Arias et al. (2021) found that the coast of Patos lagoon was the area of establishment and differentiation of main coastal lineages of the plant Calibrachoa heterophylla. This indicates that this region could have acted as a refugium for the pampean biota. The populations of Ol. nigripes haplogroup East and Ox. nasutus haplogroup North expanded inland from this centre of origin to the west and to the south. This last dispersion halted in the area between the south of Patos lagoon and the north of Mirim lagoon. In their phylogeographical study of Ox. nasutus, Pecanha et al. (2017) suggested that major elements of the Patos–Mirim lagunar complex, such as rivers and palaeochannels, could have constituted a geographical barrier to historical gene flow among haplogroups of this rodent. In line with this proposal, phylogenetic breaks for several lineages have also been reported in the area, that is for S. tumidus (Quintela et al., 2015), P. henselii (Barrasso, 2014), P. falcipes (Langone et al., 2016) and T. sidoides (Moreno et al., 2018).
The south-eastern coastal area of Uruguay was the centre of origin for Ox. nasutus haplogroup South and for S. tumidus, but at different times, 390 kya for the former and 90 kya for the latter. The relatively recent date for S. tumidus could be explained by its hydrophilic nature, which makes it prone to successive local extinctions and colonizations from few individuals. These processes, in turn, could have erased the signature of historical genetic structure, which is in line with our results for the species (Supporting Information, Fig. S4). From this area, both rodent lineages expanded following the coastline to the west and to the north until reaching the geographical barrier of the Patos–Mirim lagunar complex.
The third putative refugium is the area in the lower Paraná river, which was the centre of origin for S. aquaticus and Ol. flavescens populations, dated at 720 and 250 kya, respectively. Rosetti & Remis (2012) indicated a similar centre of origin for the eastern clade of the grasshopper Dichroplus elongatus. Near this area, Ramos-Fregonezi et al. (2017) suggested that the lower Uruguay river basin was the origin for the freshwater fish Cnesterodon decemmaculatus, which then expanded to coastal areas and Pampas water courses. Dispersal routes for S. aquaticus and Ol. flavescens from this refugium are not coincident, since the former expanded following coastal areas to the north and to the south, while the latter dispersed inland to the west, south and north, probably as a result of their different habitat preferences.
The fourth shared centre of origin is the Salado river basin in central–eastern Argentina. According to our results, this area appears as the refugium for N. lasiurus and C. musculinus populations, dated at 500 and 250 kya, respectively. For the other three refugia inferred here, some sample localities were close to the estimated centres, but for N. lasiurus and C. musculinus this is not the case: the nearest sampled locality to the proposed centre of origin in C. musculinus is 150 km to the north, and in N. lasiurus 260 km to the south. Nevertheless, it is important to note that the RRW method is based on continuous geographical coordinates and not strictly on point coordinates; thus, the estimation obtained could be the result of averaging genetic and geographical distances of the localities included in the present study according to Lemey et al. (2010). Hence, the area shown in Figure 4 represents the estimation of 80% of the highest probability density for the root location. New analyses with more exhaustive sampling could confirm this putative refugium for populations of these two species. The Salado river flows through floodplains eastwards to the Atlantic Ocean and it is possible that its extensive and flat basin was less affected by changes in sea level during the Pleistocene, resulting in a more stable environment. Furthermore, the Salado river basin has already been proposed as a centre of origin for the frog Physalaemus fenandezae (Barrasso, 2014). From the Salado river basin N. lasiurus and C. musculinus would have dispersed to the north, south and west of the Pampas in Argentina. A similar pattern was inferred in P. fenandezae (Barrasso, 2014).
CONCLUSION
In summary, the seven sigmodontine rodents from the Pampas analysed here presented discordant phylogeographical patterns and experienced asynchronous demographic processes, as postulated for several other species in the region. Four centres of origin (refugia) associated with watercourses or coastal areas were inferred: the west coastal plains of Patos lagoon, the south-east coast of Uruguay, the lower Paraná river area and an area in the Salado river basin. Most of the dispersal routes were coincident in species sharing a centre of origin, towards the north, south and west, probably following the expansion of suitable habitats in the Pampas ecoregion. The demographic stability or increase in population sizes observed in these and other pampean species from the Middle Pleistocene to the present suggest a milder impact of climatic oscillations in the region compared to those cyclically covered with ice (Patagonia) or having higher continentality (Amazonia). The ecological characteristics of the species would then have had a prominent influence on their phylogeographical patterns, producing a more complex scenario than those proposed based only on the basis of climatic oscillations.
SUPPORTING INFORMATION
Additional supporting information may be found in the online version of this article on the publisher’s website.
Figure S1. Representation of the distribution in South America (demarcated in grey; according to Patton et al., 2015) and sampling sites (represented as dots) in the Pampas ecoregion (demarcated in yellow) for each species. A, Calomys musculinus; B, Necromys lasiurus; C, Oligoryzomys flavescens; D, Oligoryzomys nigripes; E, Oxymycterus nasutus; F, Scapteromys aquaticus; G, Scapteromys tumidus. For detailed information on these species see https://cma.sarem.org.ar/, https://www.flickr.com/photos/rothfauna/9140890430/in/photostream/ and https://www.flickr.com/photos/rothfauna/9140771072).
Figure S2. Mismatch distribution analyses of each species/haplogroup. A, Calomys musculinus; B, Necromys lasiurus; C, Oligoryzomys flavescens; D, Oligoryzomys nigripes haplogroup East; E, Oxymycterus nasutus haplogroup North; F, Oxymycterus nasutus haplogroup South; G, Scapteromys aquaticus; H, Scapteromys tumidus. Grey and black lines represent simulated and observed frequencies, respectively.
Figure S3. Bayesian Skyline Plot (BSP) analyses of each species/haplogroup. A, Calomys musculinus; B, Oligoryzomys flavescens; C, Oligoryzomys nigripes haplogroup East; D, Oxymycterus nasutus haplogroup North; E, Oxymycterus nasutus haplogroup South; F, Scapteromys aquaticus; G, Scapteromys tumidus; H, Necromys lasiurus.
Figure S4. Correlation between genetic and geographical distances. A, Calomys musculinus; B, Oligoryzomys flavescens; C, Oligoryzomys nigripes haplogroup East; D, Oxymycterus nasutus haplogroup North; E, Scapteromys tumidus. Significant values are in black and non-significant values are in grey.
ACKNOWLEDGEMENTS
We are grateful to Ulyses F. J. Pardiñas, María Laura Martin, Silvana Levis and José Priotto for providing some of the specimens used in this study, and to Laura Patricia Salazar, who helped with preparation of the figures. We also acknowledge Julia Mariano, Marcelo Kittlein and the anonymous reviewers, whose comments and suggestions greatly helped to improve the manuscript. This research was supported by grants of the Agencia Nacional de Promoción Científica y Tecnológica, Argentina (PICT 2016 No. 1328 and PICT 2019 No. 01545) and of the Secretaría de Ciencia y Tecnología (Universidad Nacional de Córdoba). C.N.G., R.E.G-I. and P.C.R. are career researchers of CONICET. N.O. is a postdoctoral fellow of the Agencia Nacional de Promoción Científica y Tecnológica. The authors declare that there is no conflict of interest.
DATA AVAILABILITY
Table S1 in the Supporting Information includes the data sources used in this study.
