So Far Away, Yet So Close: Strong Genetic Structure in Homonota uruguayensis (Squamata, Phyllodactylidae), a Species with Restricted Geographic Distribution in the Brazilian and Uruguayan Pampas

The Pampas is a biologically rich South American biome, but is poorly represented in phylogeographic studies. While the Pleistocene glacial cycles may have affected the evolutionary history of species distributed in forested biomes, little is known about their effects on the habitats that remained stable through glacial cycles. The South American Pampas have been covered by grasslands during both glacial and interglacial periods and therefore represent an interesting system to test whether the genetic structure in such environments is less pronounced. In this study, we sampled Pampean populations of Homonota uruguayensis from Southern Brazil and Uruguay to assess the tempo and mode of population divergence, using both morphological measurements and molecular markers. Our results indicate that, in spite of its narrow geographic distribution, populations of H. uruguayensis show high levels of genetic structure. We found four major well-supported mtDNA clades with strong geographic associations. Estimates of their divergence times fell between 3.16 and 1.82 million years before the present. Populations from the central portion of the species distribution, on the border between Uruguay and Brazil, have high genetic diversity and may have undergone a population expansion approximately 250,000 years before the present. The high degree of genetic structure is reflected in the analyses of morphological characters, and most individuals could be correctly assigned to their parental population based on morphology alone. Finally, we discuss the biogeographic and conservation implications of these findings.


Introduction
In spite of its rich biological diversity, South America is one of the least studied continents in terms of its phylogeography [1]. Phylogeography is helpful for revealing the patterns of is an interesting model system to test these hypotheses as the strong ecological requirements of this species may limit dispersal among populations and lead to population differences in the absence of physical barriers or climatic oscillations. In this study, we take a phylogeographic approach to better understand the evolutionary history of H. uruguayensis, using molecular markers and morphological measures. More specifically, we aim to answer the following questions: How strong is the genetic structure in this species and what is the temporal depth of the mtDNA genealogy? Is there evidence of changes to this species' population size during the Pleistocene? Is morphological variation within this species concordant with the levels of population structure suggested by genetic data?  Table 1. Major mtDNA clades (see text for details) are represented by stars (Clade I), triangles (Clade II), squares (Clade III), and circles (Clade IV) and delimited geographically by a solid line coloured according to each major clade (Clade I-blue; Clade II-red; Clade III-orange; Clade IV-green). The dotted line represents the Haedo Range, in Uruguay, and "Coxilha Grande", in Brazil. Elevation is shown in grayscale according to the legend. The Uruguay River form

Sample Collection
We collected 106 individuals of H. uruguayensis from 12 localities (Table 1; Fig. 1) in Brazil and Uruguay covering the entire known distribution for this species. Geographic distances among sample sites range from 25 km to 275 km. Dr. Mariana Morando (CENPAT-Puerto Madryn, Argentina) kindly provided DNA samples from one individual of H. fasciata Duméril & Bibron 1836 and H. borellii Peracca 1897, collected from the Argentinean provinces of Santiago del Estero and Mendoza, respectively. These two species were used as outgroups.

Ethics Statement
Collection of H. uruguayensis specimens was approved by the governmental authorities of both countries (Ministério do Meio Ambiente, Brazil-SISBIO 12613-1 and Dirección de Recursos Naturales del Ministerio de Ganadería, Agricultura y Pesca, Uruguay). Euthanasia of all specimens was performed by injecting 2mg of 10% ketamine hydrochloride close to the animal's heart (http://www.asih.org/sites/default/files/documents/resources/guidelinesherpsresearch2004.pdf). The abovementioned permits include the method of euthanasia employed in this study. Specimens collected in Brazil were deposited in the Herpetology Laboratory's scientific collection at the Federal University of Rio Grande do Sul. Specimens collected in Uruguay were deposited in the scientific reptile collection of the Faculty of Sciences, at the National Museum of Natural History in Montevideo. Currently, H. uruguayensis is not designated as an endangered species in any list.

Molecular Methods
Tissue from the liver, tail muscle or tongue was used for DNA extraction. Immediately after collection, individuals were preserved in absolute ethanol and then stored at-20°C. DNA was extracted using a salting-out method as in Medrano et al. [27]. We studied two fragments of the mitochondrial DNA (720bp of the cytochrome b gene-cytb, and 420bp of the 12S rDNA gene-12S), which were amplified by Polymerase Chain Reaction (PCR) using the primers 703Botp (5'TCA AAY ATC TCA ACC TGA TGA AAY TTY GG3') and MVZ16p (5'GGC  AAA TAG GAA GTA TCA YTC TGG YTT3') as modified from Pook, Wüster & Thorpe [28], and L1091 (5'AAA CTR GGA TTA GAT ACC CYA CTAT3') and H1478 (5'GAG GGT GAC GGG CGG TGT GT3') as modified from Kocher et al. [29] for cytb and 12S, respectively. For all reactions, we used 1.5 mM MgCl2, 1X PCR Buffer, 0.2 μM of each dNTP, 0.2 μM of each primer, and 0.2 U Taq DNA polymerase. We used the cycling conditions suggested by Carranza et al. [30] with an annealing temperature of 56°C. PCR products were checked in 1% agarose gel stained with GelRed TM , purified with exonuclease I and shrimp alkaline phosphatase (ExoSAP, GE Healthcare) and Sanger sequenced in both directions at Macrogen (Seoul, South Korea).

Sequence Analysis
Chromatograms were visually inspected in the software Chromas v.2.33 (http://www. technelysium.com.au) and manually edited using MEGA 5 [31]. The cytb gene was aligned using the MUSCLE algorithm in MEGA 5 with the default parameters. The 12S gene was aligned using the strategy Q_INS-I implemented in MAFFT v6 [32], which takes into account 12S secondary structure [33]. The two mtDNA genes were concatenated in BioEdit 7.0.5.3 [34] and analysed together. All sequences generated during this study are available from GenBank (KM677689-KM677904) (S2 Table).

Phylogenetic Analysis and Divergence Times
We checked for substitution saturation using the I SSC statistic in the program DAMBE 5.2.27 [35], [36]. We did not find evidence of substitution saturation in 3 rd codon position (I SSC < I SSC, C ), and therefore all positions were used for phylogenetic analysis and estimates of divergence time. The best evolutionary nucleotide model for maximum likelihood (ML) and Bayesian inference was selected using Partition Finder v.1.1.1 [37]. Phylogenetic trees were inferred for the different haplotypes using both ML and Bayesian inference. We conducted the ML analysis in the program RAxML 7.7.1 [38] on the web server http://embnet.vital-it.ch/ raxml-bb using the GTR+G model for each of the three different partitions identified by Partition Finder: 12S gene plus 1 st codon position of the cytb gene, the 2 nd codon position of the cytb gene, and the 3 rd codon position of the cytb gene. RAxML uses GTR+G as the default model upon which simpler models can be derived depending on the parameter values that best fit the observed data. Searches were based on 100 rapid bootstrap replicates (RBS), followed by a thorough ML search. RAxML analyses were run in three replicates to check for convergence in tree topology and likelihood score. Bayesian inference was performed using the BEAST 1.8.0 software [39] with nucleotide partitions and models suggested by Partition Finder (HKY+G for 12S+cytb1 st , HKY+I for cytb2 nd , and TrN93+G for cytb3 rd ). Analyses based on different haplotypes used a Yule tree-prior and 20 million steps of the Markov Chain Monte Carlo (MCMC) sampling every 1,000 steps. This program was also used to compute the Bayes Factor to compare alternative mtDNA haplotype tree topologies based on a Path Sampling strategy [40]. We assessed chain convergence by comparing the results of at least three independent runs, and we considered MCMC sampling to be sufficient whenever ESS values were above 200 for all parameters. Convergence and the effective sample size of all MCMC runs were checked in Tracer 1.6 (http://tree.bio.ed.ac.uk/software/tracer/). The support values for the clades obtained in the Bayesian analysis are given in terms of their posterior probability (PP), which is conditioned on the trees sampled along the MCMC run. There are no fossils or known geological events for this group that could be used for a specific estimate of the substitution rate. Therefore, we used the evolutionary rate estimated by Arnold et al. [41], who analysed the same mtDNA genes for Hemidactylus, a gecko from the Gekkonidae family, and suggested an evolutionary rate of 1.15% per million years (Myr -1 ). This rate is consistent with the widely used "general reptile mtDNA rate" that was reported by Zamudio & Greene [42], which ranges from 0.47% Myr -1 to 1.32% Myr -1 . It is also consistent with estimates for other mtDNA genes reported for Crenadactylus ocellatus, a gecko from the Diplodactylidae family [43]. Divergence times were estimated in BEAST assuming the strict clock model and a normal distribution for the substitution rate, with a mean of 1.15% Myr -1 , as explained above, and a standard deviation of 0.20% Myr -1 to allow for some uncertainty in the evolutionary rate. The strict clock model is usually well justified when most data are intraspecific or from closely related species [44]. A recent phylogenetic study in Homonota failed to reject a strict clock model even with interspecific data and older divergence events compared to our study [19].

Genetic and Geographic Structure
Summary statistics, such as the number of different haplotypes, haplotype and nucleotide diversity, were estimated in the program Arlequin 3.5 [45]. This program was also used to assess the level of genetic structure among subpopulations using F ST , which is analogous to Wright's F-statistics but takes into account the genetic distance among haplotypes [46]. The evolutionary relationship among haplotypes was inferred using a median-joining network [47] in the program Network 4.6 (www.fluxus-engeneering.org). The correlation between genetic and geographic distances was tested using a Mantel test [48] in the program Alleles in Space 1.0 [49]. We tested for deviations from the null hypothesis of a constant population size using the neutrality tests of Tajima's D [50] and Fu's F S [51] in the program Arlequin 3.5 [45]. We also performed a Bayesian Skyline plot (BSP) analysis, which does not assume a priori any growth model and infers effective population size through time based on coalescent theory [52]. The BSP was estimated in the program BEAST 1.8.0 [39] as described above, except that the coalescent BSP tree-prior was used instead of the Yule model used in the previous analysis. A joint estimate of the effective population size and exponential growth rate for each population was also performed using the Bayesian search strategy in the program LAMARC v. 2.1.6 [53]. The estimate was based on three replicates of 20 short initial chains of 10,000 steps and two long chains of 1 million steps sampled every 100 steps. Initially, we estimated effective population size and population growth simultaneously to test for significant growth. As a next step, for all populations with no evidence of population growth, we estimated the effective population size based on a constant population size model using the search strategy detailed above. This strategy aimed to reduce the uncertainty in population size due to large but uninformative values for the growth parameter. We assumed the same partitions for demographic estimates in BEAST and LAMARC as in the phylogenetic analyses.

Morphological Analysis
To compare the level of genetic structure and the degree of morphological differentiation, we measured 16 traits in 79 adult individuals (35 male, 44 female) for which we had sequence data. Individuals were considered adults if they were larger than 38.47 mm for males or 35.08 mm for females (L. Martins, personal communication). We measured the following traits using a digital pachymeter with 0.01 mm precision: Snout-vent length (SVL), head length (HL), head width (HW), neck length (NL), neck width (NW), snout length (SL), base tail width (BTW), interlimb length (ILL), humerus length (HUL), forearm length (FAL), tibia length (TL), femur length (FL), length of the third and fourth fingers (3FL and 4FL), and length of third and fourth toes (3TL and 4TL).
Initially, we extracted the natural logarithm of all values and tested the interdependence of traits HL vs. HW, and NL vs. NW using a standard regression. We tested for sexual dimorphism using Student's t or Mann-Whitney test depending on the normality of data, which was tested using the Shapiro-Wilk W test. Differences in SVL among sampling sites were tested using one-way ANOVA and the post-hoc Tukey test. We used a forward stepwise discriminant analysis to estimate which variables were the most effective in assigning an individual to its original population. All analyses of morphological data were performed in the programs PAST v.2.02 [54] and STATISTICA v.7.1 (StatSoft Inc).

Results
A total alignment of 1,088 bp, including the two mtDNA genes (694 bp for cytb; 394 bp for 12S), was obtained for all 108 individuals. Overall, we identified 127 variable sites, of which only 15 were singletons, and 50 different haplotypes for the concatenated mtDNA genes in H. uruguayensis. The alignment of the 12S gene also showed five polymorphic indels in H. uruguayensis, which were treated as missing data in further analyses. Of these, two may represent synapomorphies (one deletion shared by Clades III and IV, one insertion shared by two haplotypes from Clade IV), and the other three are most likely homoplasies occurring in different branches of the mtDNA haplotype tree. We did not find substitutions causing insertion-deletion polymorphisms or stop codons for the cytb sequences, nor double peaks in the chromatograms that suggested a heterozygous individual. These results suggest that our sequences are genuine mtDNA sequences rather than inserts of mtDNA in the nucleus (NUMTs).

Phylogenetic Analysis and Time Estimates
ML and Bayesian methods resulted in different trees. In the ML tree (S1 Fig.), the H. uruguayensis clade receives only moderate support (Bootstrap = 84), and haplotypes from ART1 (all populations abbreviations follow those in Table 1) appear paraphyletic to all others. On the other hand, the Bayesian tree (Fig. 2) recovered the H. uruguayensis clade with maximum support (PP = 1.00) with a clade containing the haplotypes from the CTI sister to all others (PP = 0.99). The long branch between H. uruguayensis and the closest outgroup (H. borellii) could suggest an unstable root in the ML analysis. To test whether both topologies were equally likely, we ran a Bayesian analysis enforcing the topology found by the ML search and then compared the analyses using Bayes Factors, which indicated that the unconstrained topology found in our original Bayesian analysis is~10 37 times better supported than the alternative ML topology, which is highly significant [55]. Therefore, we based all of our data interpretation on the Bayesian phylogenetic tree.
There are four major mtDNA lineages in H. uruguayensis (which we call Clade I-IV) that diverged more than 1 Myr ago, and several well supported clades within these four major Clades (which we call "subclades") ( Fig. 2). Clade I is exclusive of the CTI population, Clade II contains the haplotypes for the Uruguayan populations RIV, PAY, and SAL, Clade III contains haplotypes found in FSF and FCB, and Clade IV contains haplotypes found in all remaining populations that are located more "centrally" in the species distribution (see Fig. 1). Only a few populations show evidence of having some phylogenetic structure. For example, BR293 had haplotypes from two different subclades; ART1 was paraphyletic in relation to ART2, and CVE also exhibited some internal structure, though not as much as BR293. Interestingly, all populations with higher phylogenetic structure are located within the central region of the species distribution.

Genetic Structure and Population History
In general, haplotype diversity was high, with an average of 0.97 for the whole sample, but there were considerable differences among populations (Table 1). We excluded FSF from all population and morphological analyses because we only sampled two individuals from that population. Populations from the central region usually had higher values (CJA, BR293, CVE, ART1), as did the divergent CTI population. On the other hand, populations such as SAL, FBC and ART2 (the latter also from the "central region") had lower values. The haplotype network (Fig. 3) shows no haplotype sharing among localities, and that 64% of all haplotypes were singletons. This finding reinforced the high haplotype diversity for this species as a whole and the deep divergence among clades. For example, there are at least 40 mutational steps between haplotypes found in CTI (Clade I) to haplotypes from the central region (Clade IV), 25 steps between Clades II and IV and 18 steps between Clades III and IV.
Such strong genetic structure is also clear from the F ST values (Table 2), which ranged from 0.24 (CVE vs. BR293) to 0.99 (SAL vs. CTI). Additionally, as expected from the haplotype network, the F ST between populations from the central region were generally lower than any other pairwise comparison, including populations from the same clade. Uncorrected pairwise genetic distances ranged from 0.8% (CJA vs. BR293) to 5.4% (SAL vs. CTI), whereas intrapopulational distances ranged from 0.0% (SAL) to 0.7% (BR293) ( Table 2). At the mtDNA clade level, the mean distance among clades was 3.7%, ranging from 2.4% (Clade II vs. Clade IV) to 4.8% (Clade I vs. Clade II), whereas distances within clades ranged from 0.3% (Clade I) to 1.4% (Clade II). The Mantel test showed that geographic distance was significantly associated with genetic distance (P = 0.0009) and accounted for approximately 40% of the variance in genetic distances (r 2 = 0.394).
At the population level, demographic estimates from LAMARC showed evidence of population growth only for CVE; at the clade level, there was evidence of population growth only for Clade IV. There was no evidence of population growth for the whole species or for Clades I-III combined ("Peripheral populations"; Table 3), suggesting that population growth in Clade IV is not an artefact of analysing more samples. It is also noteworthy that some populations from Clade IV (CVE, BR290, ART2) produced significant results in one of the neutrality tests (Table 1). In agreement with these results, the BSP for the entire species showed a pattern of  Prior bounds for population growth between-5,000 and 15,000. 2 Estimates for female effective population size assuming exponential population growth. 3 Estimates for female effective population size assuming constant population size. 4 "Peripheral" populations include Clades I, II and III. NS-Not significant. *See text for details. doi:10.1371/journal.pone.0118162.t003 Phylogeography of Homonota uruguayensis constant population size (Fig. 4A). Comparing the BSP dynamics for the Peripheral populations and Clade IV, we found no signal of a population expansion in the former population (Fig. 4B). However, the data suggest that Clade IV underwent a population expansion beginning approximately 250,000 years ago, even though credible intervals are broad (Fig. 4C). The signals of a likely population reduction in recent times are most likely artefacts of using pooled samples from structured populations, which mimics the coalescence behaviour of lineages from a shrinking population [56]. Taking into account the estimate of five years as the generation time in H. darwinii [57], the effective female population size of H. uruguayensis is approximately 350,000 individuals (

Morphological Analysis
All of the raw morphological measures are available in S1  Phylogeography of Homonota uruguayensis for males (Wilks' Lambda = 0.0018, F(108,143) = 1.8188, P = 0.0004), for which the most important traits were HL, FL, NL, and 4TL. Overall, in spite of the low sample size for each population, 82% of both males and females were correctly assigned to their original population, and populations CTI, CJA and ART2 had 100% correct assignments for both sexes. These results indicate that there are morphological differences among populations in this species.

Genetic Structure and Polymorphism
Our data reveal a strong genetic structure in H. uruguayensis, for which mtDNA haplotypes show a deep evolutionary history. The genetic distances between H. uruguayensis clades are comparable to those found between different species of the Lacertidae family, which has 2.5 to 6.4% interspecific distances for the mtDNA 12S gene [58]. Even when the comparison involves other geckos (Phyllodactylidae [59] and Gekkonidae [60]- [62]), the values observed in this study are still high. In the genus Tarentola, values for the Kimura 2-parameter distance between subspecies and species ranged from 3 to 9% for the same genes used in our study [30]. Applying this genetic distance to our data produces a 6.5% divergence between Clade I and Clade II and an approximate 5% divergence between Clade I and either Clade III or IV. Deep divergences have also been reported for some phyllodactylid geckos. For example, there are at least eight deep evolutionary lineages that may represent different species in the Phyllopezus pollicaris complex [63]. However, P. pollicaris is distributed across the whole system of open South American biomes known as the "dry diagonal" that spans from northeastern Brazil to northeastern Argentina, which complicates the comparison with a species of restricted distribution like H. uruguayensis. Nonetheless, the results found in our study, P. pollicaris and other non-Phyllodactylidae South American gecko complexes [26] indicate that a strong mtDNA structure may be a common pattern for these organisms. The strong genetic structure found for mtDNA (Fig. 2, Table 2) might reflect a difference in female/male dispersal rather than a general trend for this species. It is known that several species exhibit male-biased dispersal, which would result in high values of mtDNA structure but lower genome-wide levels (e.g., [64]). However, previous field work on this species suggests that this is not the case. Vieira [25] tracked individuals from the FCB population and found that both home range and movement patterns were higher for females than males, both in the breeding and non-breeding seasons. In any case, dispersal was very limited, as the maximum dispersal recorded in the non-breeding season was 182.52 m for a female and 81.83 m for a male. In the breeding season, these values were 190.38 m for a female and 33.61 m for a male. Such strong philopatry for both sexes is most likely associated with the high habitat specificity of H. uruguayensis, which is found exclusively in basalt outcrops used for shelter and nesting [18].

Phylogeographic Hypotheses and Evolutionary History
The evolutionary origin and initial diversification of Homonota dates back to the Early Miocene, based on an estimated 20 to 15 My origin of Homonota [19], [65]. This time period coincides with the initial phase of the rise of the Andes [66], which affected many groups of South American animals (e.g., [67]). Another important factor affecting southern South America were the three marine transgressions in the Middle to Late Miocene, which ultimately led to the formation of wide plains that were rapidly colonised by grasslands when the sea level decreased [68]. These processes were most likely important in the evolutionary history of Homonota, as they first extinguished ancestral lineages and then allowed the diversification of new lineages from ancient "refuge" areas unaffected by marine transgressions [19], [69]. Indeed, the origin of H. uruguayensis itself may be related to a marine transgression event [19].
Three major rivers in the Brazilian and Uruguayan Pampas (the Uruguay river, to the west, the Negro river to the south, and the Ibicuí river to the North) set the bounds of the species distribution, which is predicted for species with low dispersal abilities [26], [70]. These rivers may act as geographical barriers for other species. For example, it has been suggested that the Ibicuí river is associated with genetic isolation between species and populations of Ctenomys, a subterranean mammal [71], [72], and that the Negro river may act as a barrier for two subspecies of the plant Petunia axillaris [15].
Our results may be reconciled with a historical scenario in which individuals from a centrally distributed meta-population eventually disperse to colonise new areas, or a scenario in which a more highly connected population eventually loses connectivity with its peripheral populations and thus strengthens its genetic structure. We consider the second hypothesis more likely because, in the first scenario, recurrent dispersal from the same "source" population would have resulted in paraphyly of the genetic lineages in this population, as is the case for recent island species of Bothrops vipers [73]. Nevertheless, in both cases, the most important geographic feature explaining genetic and population divergence may be rivers, which may have varied in course, water flow, or riparian vegetation through time. The extension of gallery forest along the river margins appears to have changed between glacial and interglacial cycles [6].This extension might have facilitated population isolation during periods of forest expansion, at which time suitable rocky outcrops would not be available for H. uruguayensis.
In general, the distribution of mtDNA clades is concordant with the distribution of major rivers (Fig. 1). The two possible exceptions to this rule occur in Clade IV. The Ibirapuitã river, which may be a barrier separating Clades III and IV downstream, does not act as an effective barrier upstream, as populations CVE and BR293 are closely related, even though there is no haplotype sharing between them. Our results also failed to suggest that the Quaraí river is an historic barrier for H. uruguayensis, because Clade IV is distributed on both banks of this river (Fig. 1). Although the ART1 and ART2 populations (in the Uruguayan bank) form a wellsupported subclade within Clade IV, other populations from the Brazilian side show comparable values of genetic structure (Table 3). These examples may suggest that geographic distance per se may provide a more parsimonious explanation for our results, as geographic distances accounted for approximately 40% of the genetic distances in the Mantel test. Alternatively, these rivers may be better seen as an unstable barrier, allowing population connection during certain periods, but restricting gene flow on other occasions. A similar mechanism (i.e., fluctuation in flow volume) has been evoked to explain the genetic structure of lizards [74], [75] and rodents [76] across the São Francisco River, in northeastern Brazil.
The stringent ecological requirements of H. uruguayensis may explain why this species has a limited distribution despite the lack of obvious geographic barriers to the east of the species distribution. Together with the low dispersal abilities of H. uruguayensis, such stringency may have facilitated population isolation. Indeed, the basalt outcrops used for shelter and nest construction by this species may be considered ecological "islands" in the grassland "ocean" of the Pampas. Because of its low dispersal, this species is highly dependent on stepping-stones of suitable areas to connect different populations in order to maintain genetic connectivity among populations, colonise new habitats, and re-colonise areas after local extinction. Therefore, despite rivers representing reasonable barriers to gene flow, the strong population isolation within each clade may simply reflect low mobility coupled with high ecological specificity of H. uruguayensis.
Genetic diversity values and demographic estimates suggest important differences in the demographic dynamics between different regions across the species distribution (Table 1,  Table 3, Fig. 4). In general, populations from the central part of the distribution are more variable and have larger female effective population sizes, although the credible intervals are broad because only mtDNA was used in this analysis [77]. The population expansion that affected populations from the central part of the species distribution is roughly coincident with the origin of several mtDNA subclades for Clade IV (Fig. 2, Fig. 4). Given that there is no haplotype sharing and that all subclades are restricted to a single geographic population, it is possible to infer that several rocky outcrops in the central region were colonised (or re-colonised) during this expansion, with genetic drift and population isolation later restricting the number of subclades in each population. Turchetto et al. [15] also documented a population expansion for Petunia axillaris parodii that is sympatric to H. uruguayensis. However, the population expansion for P. a. parodii is more recent, around 0.1 Myr and therefore not driven by the same environmental conditions. To our knowledge, all other studies of Pampean species tend to suggest dates more recent than 0.1 Myr, either for population expansion or divergence among genetic lineages [13]- [16]. Thus, it is striking that we have estimated a population expansion in a single clade older than 0.1 Myr. These results highlight that even in a biome mostly covered by grasslands during the entire Pleistocene, demographic patterns exhibited by different species may be largely discordant.
On the other hand, common evolutionary patterns among different species may also exist. Speranza et al. [14] found high chloroplast DNA diversity for populations of the plant Turnera sidoides occurring in the Haedo Range, which is a range of hills less than 500 m in height located in the north-northwest Uruguay. These authors also suggest that this area could have served as a refuge for T. sidoides, even though they did not analyse plants from Brazil, where the species also occurs. For H. uruguayensis, two populations with high nucleotide diversity (Table 1) and effective female population sizes (Table 3) were CVE and BR293, which are located close to the "Coxilha Grande", the continuation of the Haedo Range in Brazil. In addition, RIV in the Uruguayan portion of the Haedo Range also has high genetic diversity. This observation may support the role of the Haedo Range as a region of high genetic diversity for several Pampean species.

Morphology and Conservation
There are also morphological differences among populations, with 82% of all individuals being correctly assigned to their original population based on morphological data alone. For some populations, all individuals were correctly assigned. As the ecological context for all populations is similar, the strong population structure for H. uruguayensis may suggest that the morphological differences among populations are due to genetic drift [78]. However, it is controversial that genetic drift may be strong enough to result in morphological change in natural populations [79]- [83], and therefore we cannot exclude the role of local adaptive processes affecting H. uruguayensis morphology. It is possible that the CTI population, which has the most divergent mtDNA clade and 100% correct morphological identifications, reflects an evolutionary lineage distinct enough to be considered a different species or subspecies. However, a more comprehensive morphological characterisation and the study of variable nuclear genetic loci are necessary to answer this question. Finally, another implication of the deep population structure is that each local population contains unique genetic variation, and therefore can be considered at least as independent management units (MU) sensu Moritz [84], which is important for conservation plans for H. uruguayensis.