Cryptic lineages and Pleistocene population expansion in a Brazilian Cerrado frog

Diversification of South American species endemic to open habitats has been attributed to both Tertiary events and Pleistocene climatic fluctuations. Nonetheless, phylogeographical studies of taxa in these regions are few, precluding generalizations about the timing and processes leading to differentiation and speciation. We inferred population structure of Hypsiboas albopunctatus, a frog widely distributed in the Brazilian Cerrado. Three geographically distinct lineages were recovered in our phylogeny. The Chapada dos Guimarães (CG) clade was the first to diverge from other populations and contains multiple haplotypes from a single population in western Cerrado, probably representing a cryptic species. The southeast clade (SE) includes populations along the southeastern limit of the range within the historical distribution of the Brazilian Atlantic forest. Finally, the Central Cerrado (CC) group includes haplotypes from the interior of Brazil that are paraphyletic relative to the SE clade. Analyses of historical demography indicate significant population expansion in the CC and SE populations, likely associated with colonization of newly formed open habitats. The divergence of populations in the CG clade occurred in the late Miocene, concordant with the uplift of the central Brazilian plateau. Divergence of the SE clade from the CC occurred during the mid‐Pleistocene. Thus, both Tertiary geological events and Pleistocene climatic fluctuations promoted divergences among lineages. Our study reveals a complex history of diversification in the Cerrado, a morphoclimatic domain highly threatened because of anthropogenic habitat alteration. We identified surprisingly deep divergences in a widely distributed frog, indicating that the Cerrado is not a barrier‐free habitat and that its diversity is likely underestimated.


Introduction
In the last 20 years, advances in statistical phylogeography, focusing on the spatial and temporal changes in genetic variation within and among closely related species, have provided unprecedented details of the historical processes underlying diversification (Avise 2009;Knowles 2009). Individual species-specific studies pro-vide information on how geological, geographical, and environmental changes have altered the demographic history of populations. Combining those demographic inferences across multiple species (Arbogast & Kenagy 2001) or with data on spatial changes in habitat distribution or known barriers to gene flow (Knowles 2001;Carstens & Knowles 2006) allows us to test increasingly refined hypotheses about the history of particular species groups (Knowles 2009). In this study, we infer the evolutionary history of a Neotropical frog widely distributed in the Brazilian Cerrado and provide a microevolutionary perspective on the order and timing of divergences in this species complex. The highest global diversity of amphibians is found in the Neotropics, a pattern that has been attributed to the large number of morphoclimatic domains present in Middle and South America (Duellman 1999) each with a highly endemic and locally adapted fauna. Nonetheless, detailed knowledge of the population genetic processes leading to rangewide diversification are limited to few amphibian species that occur primarily in forested biomes (Carnaval et al. 2009;Fitzpatrick et al. 2009;Brunes et al. 2010;Thomé et al. 2010).
The main open formations in South America include the Llanos in Venezuela, the Chaco in Paraguay and Bolivia, and the savannas of the Brazilian Cerrado (Silva 1997). The Cerrado is the second largest South American domain, with a total historical area of approximately 2.5 million square kilometres (Silva et al. 2006). Most of the region is composed of seasonal savannas with corridors of mesic gallery forests. This open habitat in central Brazil is surrounded by other important forested and non-forested formations: the plains of the Amazon basin to the north, the Caatinga to the northeast, the Chaco-Pantanal to the southwest, and the montane Atlantic forest to the southeast (Silva & Bates 2002), therefore serving as an important link among them. The Cerrado is the largest, most species-rich and probably most threatened tropical savanna in the world (Silva & Bates 2002), and is included among the 25 most threatened hotspots on Earth (Myers et al. 2000). This open formation has been intensively modified in the last 40 years, primarily because of urban and agricultural development; currently, only 20% of its original area remains and approximately 6% of that area is legally protected (Myers et al. 2000).
The Cerrado (sensu Silva & Bates 2002) has relatively low topographic complexity, with vast plateaus separated by peripheral depressions (Silva 1997). The plateau surfaces are composed of rocks that were moulded during cycles of erosion in the Cretaceous, Mid-Tertiary, and Late Tertiary (Silva 1997;Colli 2005 and references therein). From the Late Tertiary to the Early Quaternary, these plateaus were uplifted to their present-day elevations (500-1700 m; Silva 1997;Nores 2004;Rossetti et al. 2005), causing the formation of depressions between plateaus and fragmenting the landscape. During most of the Quaternary, these depressions expanded progressively, as plateau surfaces and escarpments became eroded (Cole 1986;Colli 2005 and therein). In addition, dry and moist periods alternated from late Tertiary throughout the Pleistocene, leading to contractions and expansions of forests and open habitats on this landscape (Webb 1978;Behling 1998). Thus, the current Neotropical savannas in south-eastern and mid-western Brazil are relicts of a larger continuous area that existed during the dry conditions of the last Pleistocene glacial period (Behling & Hooghiemstra 2001). Given their potentially long history, and the dynamic changes in landscape and habitat distribution, species of open habitats have the potential to elucidate the history of the open habitats themselves, as well as that of the adjacent forests (Wü ster et al. 2005;Moraes et al. 2009;Werneck 2011).
The origins of the Cerrado herpetofauna are still poorly understood. Fossils of amphibians and reptiles from open habitats of central South America are rare, and most evolutionary studies are based on fossils from southern parts of the continent (reviewed in Colli 2005). Diversification of open habitat vertebrates in South America began during the Tertiary (Colli 2005;Garda & Cannatella 2007) when open habitats and dry climates dominated the continent (Pascual & Jaureguizar 1990;Hewitt 2004). Despite the diversity and potentially deep history of open-habitat taxa, molecular studies of divergence and speciation in the Neotropics have focused mainly on the Amazon Basin and the Atlantic forest (Lessa et al. 2003;Grazziotin et al. 2006;Carnaval & Bates 2007;Brunes et al. 2010) with a bias towards studies of mammals and birds (Salazar-Bravo et al. 2001;Nahum et al. 2003;Eberhard & Bermingham 2004). One reason for these biases is the common assumption that savannas are less diverse because of lower habitat heterogeneity (Lamotte 1983). Nonetheless, some hypotheses about regional diversification have been proposed (Werneck 2011;Werneck et al. 2011). For example, geologic and climatic events during the upper Tertiary, including the uplift of the central Brazilian plateau, are invoked to explain the vicariant distribution of South American vertebrates (Werneck & Colli 2006;Almeida et al. 2007;Maciel et al. 2010). More recent studies on mammals, plants and invertebrates have shown that both Tertiary events and Pleistocene climatic cycles contributed to diversification of lineages in the Cerrado, as well as other non-forested formations in South America (Almeida et al. 2007;Ramos et al. 2007;Moraes et al. 2009). However, phylogeographical studies on central South American anurans are rare (Camargo et al. 2006), precluding generalizations about processes leading to high diversification and the timing of those events.
The Cerrado is home to approximately 204 anuran species, most of which are widely distributed (Valdujo 2011). Cerrado anurans are considered to be habitat generalists with a high colonization capacity (Haddad & Prado 2005). However, a recent biogeographical analysis of species distributions revealed that anurans are not randomly distributed in central Brazil, suggesting that biomes adjacent to the Cerrado may have influenced diversification in this domain (Valdujo 2011). For instance, some taxa occur primarily in northern Cerrado, bordering the Amazon basin, others are distributed in the eastern Cerrado and share faunistic similarity with the Atlantic forest. Here, we report on a population genetic study of the hylid frog Hypsiboas albopunctatus, widely distributed and typical of the Brazilian Cerrado (Caramaschi & Niemeyer 2003;Valdujo 2011). This species is common in open habitats and adjacent forests in central South America (Fig. 1), ranging from central and southeastern Brazil, to northeastern Argentina, Uruguay, eastern Bolivia, and Paraguay (Frost 2011). Although commonly found in disturbed areas and considered a habitat generalist, H. albopunctatus individuals are generally found breeding close to water bodies with some water circulation (Muniz et al. 2008) in gallery forests at mid-elevation sites (400-1000 m; C. P. A. Prado, personal observations) and are rarely found below 350 m. Because of its wide distribution in central Brazil and at the ecotone between Atlantic forest and eastern Cerrado, we chose this species as a model for investigating processes and timing of diversification in the savannas of central South America. We used mitochondrial and nuclear markers to assess the distribution of genetic variability among rangewide populations of H. albopunctatus and elucidate their evolutionary history. Specifically, we tested the following hypotheses: (i) Tertiary events (e.g. uplift of the plateaus and landscape fragmentation) and climatic fluctuations causing expansions and contractions of open habitats in the Quaternary led to lineage diversification and have affected the distribution of genetic diversity in this species; (ii) given that H. albopunctatus is primarily distributed in the Cerrado, inhabiting typical open vegetation formations, the earliest divergences of this species occurred among ancestral Cerrado populations in central Brazil, followed by expansion into southeastern Atlantic forest; and (iii) given this species' habitat preferences, populations are genetically structured according to main geological plateaus within its distribution, because of isolation among regional populations by large topographic depressions.  & Haddad 2000), and six populations were sampled from the Cerrado, in central Brazil (Silva & Bates 2002). We chose Hypsiboas multifasciatus and Hypsiboas raniceps as successively distant outgroup taxa for phylogeographical analysis (Faivovich et al. 2005). Voucher specimens have been deposited in the Coleção de Anfíbios CFBH, Universidade Estadual Paulista-UNESP, Rio Claro, São Paulo, Brazil.

Molecular methods
We extracted whole genomic DNA from ethanolpreserved liver, muscle, or toe clips using DNeasy Tissue kits (Qiagen, Valencia, CA, USA), following manufacturer's protocols. Genomic DNA was eluted in 200 lL buffer for use as template in polymerase chain reaction (PCR) amplifications. We amplified and sequenced three mitochondrial gene fragments: (i) the control region and a short segment of the adjacent cytochrome b gene (hereafter referred to as control region fragment, 1024 bp), (ii) the complete NADH dehydrogenase subunit 1 and the adjacent tRNA Leu and tRNA Ile , and partial sequences of the flanking tRNA Gln  We used the same sequencing primers used in the original amplifications, with the exception of internal sequencing primers designed for the ND1 and ND2 fragments (Table 2). We sequenced each gene in both directions to avoid base-calling ambiguities. Products were column purified to remove non-incorporated terminator dye using SephadexÔ G-50 and electrophoresed on an ABI 3100 Genetic Analyzer (Applied Biosystems). Electropherograms were checked by eye and contiguous sequences for each fragment were assembled using SEQUENCHER version 4.1 (GeneCodes, Ann Arbor, MI, USA).

Bayesian population-level phylogeny
We aligned sequences for each gene separately using CLUSTALW (Thompson et al. 1994) in the MEGALIGN version 6.1.2 program of the Lasergene sequence analysis software (DNASTAR, Inc., Madison, WI, USA). For each gene, sequences were aligned using four different gap penalties (5,10,15,20) and gap length maintained constant (0.20) to identify regions of ambiguous homology (Gatesy et al. 1993). All other parameters were set to default settings. Positions that varied in alignment across this range of parameters were excluded as characters in our phylogenetic analyses. We concatenated all three mtDNA fragments for each individual and used unique haplotypes (for all genes combined) in phylogenetic analyses. We used PAUP* version 4.0b10 (Swofford 2000) and the PAUPUP version 1.0.2 (Calendini & Martin 2005) to determine the number of unique haplotypes. Population-level phylogenetic trees Rhodopsin GTAGCGAAGAARCTTCAAMGTA Bossuyt & Milinkovitch (2000) were reconstructed using Bayesian inference (BI) implemented in MRBAYES version 3.1.2 (Huelsenbeck & Ronquist 2001). We estimated the best model of evolution for each data partition using MODELTEST version 3.7 (Posada & Crandall 1998). We then used a partitioned BI search in a Markov Chain Monte Carlo analysis, run in two parallel chains, each for 5 000 000 generations, sampled each 100 generations, and excluding the first 12 500 trees as burn-in. Adequate burn-in was determined by examining a plot of the likelihood scores of the heated chain for convergence and stationarity. We used the clades recovered in the BI topology as a guide for grouping haplotypes for subsequent analyses focused on lineage divergence and clade expansion. For estimates of diversity indices and historical demographic parameters, which require comparison of geographical, rather than phylogenetic units, we grouped collecting sites according to the morphoclimatic domain where they were located, either Atlantic forest sensu lato (populations 1-11 and 14) or Cerrado (populations 12-18; Table 1, Fig. 1).

Molecular diversity and haplotype network
We estimated haplotype diversity (h) and nucleotide diversity (p) using ARLEQUIN version 3.01 (Schneider et al. 2000) with parameters set to default settings. We estimated genetic diversity indices considering each mtDNA and nuclear gene fragment separately and for the three mtDNA fragments combined. For the nuclear genes, we resolved haplotypes from heterozygous individuals using the PHASE algorithm implemented in DNASP v. 5 (Librado & Rozas 2009) with default settings. Diversity indices were estimated at two scales: range wide (including all populations) and by geographical region.
Overall mtDNA differentiation among sampled populations was quantified using pairwise F-statistics. Empirical F ST values were compared with a null distribution of no difference between the populations to test for significance. Analyses were implemented in ARLE-QUIN version 3.01, with 10 000 permutations. We used the Tamura method for molecular distance estimates (Tamura 1992) which allows for unequal nucleotide frequencies. The transition-to-transversion ratios, as well as the overall nucleotide frequencies, were computed from the original data. We assessed how genetic difference was partitioned among populations and regions using a hierarchical analysis of molecular variance (AMOVA, Excoffier et al. 1992), with 10 000 permutations and the Tamura method for molecular distance estimates (Tamura 1992). We estimated AMOVA using the three regional haplotype clades we identified in our mtDNA phylogenetic analyses.
Finally, we constructed haplotype networks for the three concatenated mtDNA genes and the two nuclear genes independently; using the median-joining method (MJN) (Bandelt et al. 1999) implemented in the program NETWORK version 4.1.0.8 (http://www.fluxus-engineering. com) with all parameters set to default settings. Nuclear haplotype networks were constructed using the haplotypes inferred in PHASE.

Genetic signature of past demographic processes
Mismatch distributions of pairwise differences between mtDNA haplotypes provide information about historical demography of populations (Rogers & Harpending 1992). Populations with historically stable demography have multimodal mismatch distributions, whereas a unimodal distribution (approximately Poisson) suggests a recent demographic expansion, possibly following a bottleneck. We conducted a mismatch distribution test using ARLEQUIN version 3.01, and calculated the expansion parameters s, h 1 , and h 0 . The fit of mismatch distributions to the expansion model was assessed with 500-1000 bootstrap replicates. Significance was estimated by comparing the sum of square deviations (SSD) between the observed data and those from the simulated expansion model; a significant P-value rejects the fit of the data to the expansion model. We also calculated Tajima's D (Tajima 1989) as a second test of deviation from neutrality; historical population growth predicts significantly negative D values (Tajima 1989). The significance of deviation from values expected under demographic stationarity was tested with 10 000 bootstrap replicates. Because the clades recovered in our phylogenetic analysis may have different histories, we considered each separately in analyses.
We estimated historical demographic parameters under a model of isolation with migration in the program IM (Nielsen & Wakeley 2001) using the mtDNA sequences. IM uses a Bayesian coalescent framework to estimate demographic parameters for a pair of populations. Historical demographic parameters were estimated for two geographical population groups: the Central Cerrado and the Atlantic forest. We ran an IM analysis for the divergence of the southeast (SE) and CC populations, by pooling all populations in each region. We did not perform IM analyses for the divergence of the CG from CC ⁄ SE populations (node A, Fig. 2) because those two groups are anciently diverged and have reached reciprocal monophyly, thus are not appropriate for the isolation with migration model. For the CC and SE populations, we estimated time since divergence (t = tl), the effective size of each population (h 1 = 4N 1 l and h 2 = 4N 2 l), the ancestral effective size (h A = 4N A l), and the migration rates between regions  (m 1 = m 1 ⁄ l and m 2 = m 2 ⁄ l). We performed a series of initial searches varying MCMC search parameters to assure adequate mixing and convergence. Final runs consisted of 10 coupled chains heated geometrically (heating parameters: g1 = 0.95 g2 = 0.99), an initial burn-in of 300 000 steps, and sampling every 10 steps for 10 000 000 steps. Each analysis was repeated three times with different random seed numbers to ensure that the process reached stationarity and converged on the same parameter estimates. We compared parameter posterior probability densities from each of the runs to determine the parameters (mutation rate corrected) of interest; posterior probabilities (HiPt) and the lower and upper bounds of the 90% highest posterior density interval (HPD90Lo and HPD90Hi) were estimated from the posterior distribution for each parameter. Comparison of the values of t and h A informs on the relative age and the ancestral population of the phylogenetic divergence event. We also compared posterior probabilities within runs for unequal effective population sizes (h 1 > h 2 ), population expansion within regions (h 1 > h A ; h 2 > h A ), and asymmetrical migration (m 1 > m 2 ) and evaluated those statistically by comparing the proportion of times one parameter was larger than the other over the course of the run. Once estimated, the mutation-scaled parameters were transformed to demographic parameters by using the arithmetic mean of mutation rates for two of the mtDNA genes included in the analyses, based on independent calibrations from studies of other anurans: 0.957% and 1.46% per lineage per million year for ND2 and control region, respectively (Macey et al. 1998;Crawford 2003).

Time to most recent common ancestor (Tmrca)
We estimated divergence times for the four major nodes in our phylogeny using a Bayesian Markov Monte Carlo method implemented in BEAST version 1.4.7 (Drummond et al. 2002;. To explore the effects of model selection on dating estimates, we implemented three coalescent models (Bayesian skyline, exponential and expansion models). Lacking fossil calibrations for any South American hylid frogs, we also compared results based on two published divergence rate estimates, the same mutation rates used in the IM analyses. We performed two independent analyses using concatenated mitochondrial sequences (control, ND1 and ND2) and assigning to the entire concatenated sequence the rates of evolution estimated for the control region and ND2 gene fragments (Macey et al. 1998;Crawford 2003). In total, we estimated the time to most recent common ancestor (Tmrca) at each clade node in six independent runs, each using a fixed mutation rate corresponding to one of the two genes for which these data were available, and one of the three possible models of evolution. For each run, we examined log files using the software package TRACER version 1.4  to determine stationarity of the posterior distributions for all model parameters and obtain average divergence time and 95% confidence intervals. Rates were uncorrelated among lineages in the phylogeny and priors were specified for the respective substitution models as follows; rate parameters uniform [0,500]; alpha exponential (1); proportion of invariant sites uniform [0,1]; default values for all scale operators. For each analysis, we performed two independent runs of 20 million generations sampling every 1000th generation and removing 10% of the initial samples as burn-in.

Bayesian population-level phylogeny
We aligned a total of 3557 bp combining the three mtDNA gene fragments for 121 individuals, after excluding 158 bp because of ambiguous homology. We identified 91 unique haplotypes across the 18 sampled localities. Duplicate haplotypes were always from individuals collected at the same site (Table 1). We identified the best-fit model for all mtDNA coding genes (cytb, 16S, ND1 and ND2), the control region, and the seven complete or partial tRNAs included in our data set. For ND2 and 16S, the best fit models were the General Time Reversible model (GTR) including invariable sites and rate variation among sites (GTR + I + G) and only invariable sites (GTR + I), respectively. For ND1, cytb, and control region, the best model was the Hasegawa model (HKY), including invariable sites and rate variation among sites (HKY + I + G) for control region and cytb, and only rate variation among sites (HKY + G) for ND1. The tRNA sequences best fit the K80 model (Leu), the HKY Model (Gln and Iln) and the HKY + I model (Trp, Ala, Asn, Cys, and Tyr). The fully partitioned Bayesian analysis of unique haplotypes for all mtDNA genes resulted in a 50% majority rule consensus tree (Fig. 2). The most obvious feature of our topology is the well-supported division of haplotypes into three widely distributed groups that are generally concordant with geography. The first lineage includes all haplotypes from the Chapada dos Guimarães (hereafter the CG clade), the western most population we sampled in the state of Mato Grosso (Fig. 1), which diverged early in the history of this species group (node A, Fig. 2) and is the sister clade of the remaining lineages. The second well-defined clade includes populations distributed along the southeastern limit of the range (hereafter the SE clade) in open habitats inside and at the edge of the Atlantic forest and transitional areas, located primarily in the states of Minas Gerais, São Paulo, Paraná, and Santa Catarina (node D, Fig. 2). The third is the Central Cerrado, which is not a welldefined clade (hereafter the CC haplotypes), but rather a group of haplotypes and subclades that form a polytomy from within which the SE clade is derived (node C, Fig. 2). These haplotypes occupy the Cerrado in the interior of the country, in the states of Goiás and Mato Grosso do Sul, with a few haplotypes from the southeastern states of Minas Gerais and São Paulo.
The three haplotype groups differ in degree of resolution and nodal support (Fig. 2). The well-supported CG clade is the most isolated, and haplotypes from this single population are highly divergent from each other (node B, Fig. 2). The SE clade is also well supported (96% posterior probability) and contains a single sister haplotype from Três Lagoas, Mato Grosso do Sul (h86). Other haplotypes sampled from this same population are included in the CC group of haplotypes (Fig. 2, Table 1). Moreover, the SE clade contains three subclades: the first includes samples from central and western São Paulo and from the highlands of Serra do Cipó , Minas Gerais state (SE1), the second subclade is composed of haplotypes primarily from sites located at the southern Atlantic forest mountain ranges, including southern São Paulo, Paraná, and Santa Catarina states (SE2), and finally, the third is composed of samples from northeastern São Paulo and southern Minas Gerais states (SE3), along the Atlantic forest mountain range of the Serra da Mantiqueira. In contrast, the CC haplotypes show low levels of population genetic structure, as evidenced by the shorter branches in the BI phylogram, lower posterior probabilities, and a basal polytomy. Nonetheless, it is possible to recognize two regional haplotype groups: the first composed of haplotypes from the central states of Mato Grosso do Sul and Goiás (CC1) and the second including a group of haplotypes from the southeastern populations of Cristina ⁄ Itamonte (h66, h68, h70, h67, h70 and h76) and Serra do Cipó (h88), Minas Gerais state, and a single haplotype from Rio Claro (h12), São Paulo state (CC2). Most haplotypes from these southeastern populations (SP, MG), however, appear in the SE clade (Fig. 2, Table 1). In summary, in a few cases, we found haplotypes that fell outside the geographical distribution of clades. These cases only occurred between southeast and central Cerrado populations; no haplotypes were shared between the Chapada dos Guimarães and other localities.

Molecular diversity and haplotype network
Based on the mtDNA sequences, the average haplotype diversity was high both within regional population groups and among populations, ranging from 0.96 in the CG clade to 0.99 across all populations (Table 3). Nucleotide diversity was also high, ranging from 0.014 in the CC clade to 0.044 in the CG clade. The high nucleotide diversity found in the Chapada dos Guimarães population was unexpected, as our sample included only 11 individuals collected from a single site. The molecular diversity indices were consistent for all three mtDNA fragments analysed (Table 3). Considering all populations, the fragment ND1 was less diverse compared with the others, whereas the control region exhibited the highest diversity, number of polymorphic sites, and mean number of pairwise differences. In contrast, the two nuclear gene fragments showed low polymorphism among populations. We resolved 12 unique haplotypes among 38 individuals for b-fibrinogen (674 bp; N polymorphic sites = 9; p = 0.001 ± 0.001; h = 0.40 ± 0.10) and 11 among 38 individuals for rhodopsin (358 bp; N polymorphic sites = 7; p = 0.003 ± 0.002; h = 0.74 ± 0.04).
Pairwise population estimates of F ST based on the mtDNA data ranged from 0.13 (Rio Claro, SP -Ribeirão Branco, SP) to 0.89 (Pindamonhangaba, SP -Piraí do Sul, PR) (F ST values are available as Supporting information). Most (114 of 153) F ST values were significant (a £ 0.05), indicating high inter-population genetic diversity. Non-significant pairwise comparisons were primarily observed between geographically close populations. Because of the high genetic divergence of the CG clade, we carried out two AMOVA analyses: (i) including the three phylogenetic groups (SE, CC and CG) and (ii) including only the two most recent phylogenetic groups (SE and CC). The two AMOVAs showed concordant results and underscore the depth of genetic divergence within and among H. albopunctatus populations. Genetic differences among regional population groups explained most of the total variance (53.5-73.5%, Table 4) independent of whether the CG clade was included or not.
The median-joining network of the concatenated mtDNA genes underscores the nonrandom association of haplotypes with geography and is fully concordant with patterns obtained in our phylogenetic analyses (Fig. 3). Our network supports the distinctiveness of the population from Chapada dos Guimarães, which was separated from the others by 187 mutational steps, whereas SE and CC haplotypes differ by 22 mutational steps. A group of haplotypes collected in southeastern populations is divergent but associated with the CC; these are the same haplotypes included in the sub-group of haplotypes CC2 in our tree (h12, h66, h67, h68, h70, h73, h76 and h88). The same single CC haplotype (h86), from Três Lagoas, Mato Grosso do Sul state, was included in the southeastern network group (Fig. 3).
We constructed median-joining networks separately for each of the two nuclear fragments (b-fibrinogen and rhodopsin). Nuclear markers were less informative than the mtDNA haplotypes, and unique alleles were generally not concordant with geography (Fig. 4). We recovered 12 unique b-fibrinogen alleles; 39 of the 76 alleles sequenced (38 diploid individuals) were identical, and that commonest allele was found in all three regions (Fig. 4). Among rhodopsin sequences, we detected 11 unique haplotypes among the 76 sequences; seven of those were shared among SE and CC regions. In both nuclear markers, alleles were typically separated from other alleles by only one mutational step.

Genetic signature of past demographic processes
Mismatch analyses significantly supported the hypothesis of range-wide population expansion, but patterns differed among the groups examined (Table 5, Fig. 5). We detected a signature of rapid population expansion both for the SE clade and CC group of haplotypes using the goodness-of-fit test; Tajima's D was significant only for the CC group (Table 5). Neither mismatch analyses nor Tajima's D identified a signature of population expansion in the CG clade.
Our IM analyses indicated that the SE and CC regional populations diverged approximately 89 000 years    Fig. 3 Medium-joining networks among mtDNA haplotypes of Hypsiboas albopunctatus populations underscore the deep genetic structure among the three geographical regions sampled. Haplotypes from southeastern Brazil populations that fall within the Central Cerrado haplotype group appear with a dark grey background. Each circle represents a unique haplotype with size proportional to its relative frequency (from 1 to 5). Bars along branches connecting haplotypes represent mutational steps. larger than the ancestral population size (N A = 70 233, CI = 41 693-117 500; proportion of time q CC > q A = 1.00, q SE > q A = 1.00; Table 6). Current population sizes for CC and SE are not significantly different (proportion q CC > q SE = 0.32; Fig. 6). Finally, we found evidence for limited, yet asymmetrical, migration between the CC and SE regions. Estimated effective migration rates (adjusted for population size) were very low (ranging from 0.11 to 2.51 frogs per generation), but historical migration into the SE populations from CC populations was significantly higher than the reverse direction (proportion of time 2N 1 m 1 > 2N 2 m 2 = 0.95).

Time to most recent common ancestor (Tmrca)
The estimates of Tmrca for the four nodes of interest in our phylogeny (A-D, Fig. 2) were highly concordant among runs assuming different models of evolution (Table 7). As expected, divergence estimates varied significantly with different gene-specific mutation rates, underscoring the importance of mutation rate in divergence estimates that are not fossil-calibrated. Even so, confidence intervals of our estimates of Tmrca for two different mtDNA genes under different models are overlapping (Table 7). Using the ND2 mutation rate and the estimates under the coalescent expansion model, the age estimate of the oldest divergence (node A) including all Hypsiboas albopunctatus haplotypes was approximately 5.8 million years ago (Ma), a date that puts the initial diversification of this group in the Miocene ⁄ Pliocene boundary. The haplotypes of the Chapada dos Guimarães clade also show relatively old coalescence dates, estimated at approximately 2.5 Ma in the late Pliocene. Finally, the populations in the widespread CC ⁄ SE clade show very recent coalescence dates (Node C; c. 0.9 Ma, in the mid-Pleistocene) and those dates are not statistically distinguishable from the divergence of the youngest clade of populations in southeastern Brazil (SE, Node D: c. 0.7 M, also in the mid-Pleistocene). The Tmrca estimates are older than the divergence time estimated for node C derived from IM, and the confidence intervals for both estimates are not overlapping (Tables 6 and 7).  Fig. 4 Medium-joining networks among nuclear (b-fibrinogen and rhodopsin) haplotypes of Hypsiboas albopunctatus populations. Both nuclear markers show low allelic diversity that is widely geographically distributed. Each circle represents a unique haplotype with size proportional to its relative frequency (indicated by numbers). Each haplotype is separated by one mutational step.

Table 5
Results of the model of population expansion using a goodness-of-fit test (Rogers & Harpending 1992;Schneider & Excoffier 1999) and Tajima's D (Tajima 1989) for Hypsiboas albopunctatus population groups. The parameters of the model of population expansion, s, h 0 and h 1 are the age of expansion and population size before and after expansion, expressed in units of mutational time. We report goodness-of-fit tests for the mismatch distributions (based on SSD, sum of squared deviations) and significance of D

Population-level phylogeny and divergence times
Our analyses revealed deep divergences among three regional lineages of Hypsiboas albopunctatus. The genetically diverse CG clade is evidence of a cryptic lineage with an unusually high degree of genetic diversity, considering that samples were collected from a single locality. Finer scale assessment of the distribution of this lineage, as well as potential contact zone between the cryptic form and populations within the CC haplotype group, will be necessary to better understand diversification in this region. The second well-supported clade in our topology is composed of haplotypes from the southeastern part of the species' range (SE) that seemingly evolved from a group of haplotypes in central Cerrado (CC). Our Bayesian divergence time estimates placed the divergence between CG and CC at 5.8 Ma, during the late Miocene (node A, Fig. 2). These estimates are highly concordant with the uplift of the central Brazilian plateau (Saadi et al. 2005), which may also have contributed to old divergences in other Neotropical vertebrates, including mammals (Almeida et al. 2007), birds (Silva 1995), reptiles (Werneck & Colli 2006) and amphibians (Garda & Cannatella 2007;Maciel et al.  for the lineage splitting event that gave rise to the CC and SE populations of Hypsiboas albopunctatus (node C, Fig. 2). Parameter estimates are scaled by mutation rate; the demographic parameters inferred from these estimates are reported in Table 6.  Table 5. 2010). Thus, our phylogenetic analysis and estimates of divergence times corroborate the hypothesis of a central Brazilian origin for H. albopunctatus with deep divergences early in the history of the species group, followed by a colonization of the adjacent Atlantic forest in southeastern and southern Brazil. The uplift of the Chapada dos Guimarães plateau occurred at the end of the Tertiary, during the Plio-Pleistocene transition (Silva 1997;Porzecanski & Cracraft 2005), coinciding with the divergence time estimated for haplotypes in the CG clade (node B, Fig. 2). During that period, dry climate persisted for millions of years (Ab'Sáber 2006). The Chapada dos Guimarães has a complex topography; the region is located at the edge of a large plateau approximately 600-800 m in elevation, with an abrupt border at the southern end, dropping abruptly to 200-300 m (Ab'Sáber 2006). In an analysis of regional differences within the Cerrado, the Chapada dos Guimarães was included in ecological unit 1B (Silva et al. 2006), a unit characterized, among other traits, by elevations between 550 to 700 m. Although this unit is patchily distributed throughout the Cerrado, the area of the Chapada dos Guimarães plateau is almost completely isolated, surrounded by lower elevations (70-400 m). Thus, populations remaining in mesic vegetation patches there were likely isolated, promoting interruption of gene flow, increase in drift, and fixation of genetic variation. During our collections, we never registered H. albopunctatus below 350 m of elevation in this region (C. P. A. Prado and K. R. Zamudio, personal observations), because at lower elevations this species seems to be replaced by its sister taxon Hypsiboas raniceps. Two other frog species, Phyllomedusa centralis and Pristimantis crepitans, are known only from type localities in the Chapada dos Guimarães and adjacent plateaus (Frost 2011), underscoring the isolation of this region relative to other mid-elevation plateaus in central South America.
The diversification of Central Cerrado haplotypes occurred approximately 0.9 Ma followed by the divergence between the CC and the SE clade 0.7 Ma (nodes C and D, respectively), both during the mid-Pleistocene. Moist and dry climate alternated between the late Pliocene throughout the Pleistocene, repeatedly changing the historical distribution of the Cerrado (Webb 1978;Behling 1998;Behling & Hooghiemstra 2001). It is likely that the CC populations and SE clade became differentiated during dry periods, when Cerrado vegetation expanded into southeastern Brazil. These estimates of Table 6 Historical demographic parameters for Hypsiboas albopunctatus lineages in the Brazilian Cerrado. We inferred historical demographic parameters for populations grouped according to geographic region (SE and CC are Populations 1 and 2, respectively). Mutation-scaled parameters derived from IM were converted to demographic parameters using the arithmetic mean of mutation rates per million years for two mitochondrial genes (ND1 and control region; see text) 2.87 · 10 )6 8.02 · 10 )7 -6.54 · 10 )6 Effective migration 1 ‡ 2.506 0.548-7.255 Migration 2 † 1.15 · 10 )7 1.15 · 10 )7 -2.18 · 10 )6 Effective migration 2 ‡ 0.112 0.083-2.785 *Number of individuals. † Number of migrants per generation. ‡ Population migration rate (mutation-scaled migration multiplied by population size), the effective rate at which genes move into populations per generation.  (Haffer 1997). During Pleistocene dry periods, gallery forests may have played an important role in diversification and persistence, by providing refuge for riparian-associated species (Prance 1973) and promoting migration and genetic contact among otherwise isolated areas. A number of plant species of the Amazon forest are distributed throughout the gallery forests in central and southeastern Brazil, corroborating the importance of this habitat for the persistence of some taxa (Prance 1973). Hypsiboas albopunctatus is adapted to open habitats, but has a preferred elevational range (c. 400-1000 m; Fig. 1) and is generally associated with gallery forests and small streams or ponds with some water circulation (Muniz et al. 2008). If current habitat use by this frog is any indication of ancestral preferences, then gallery forests may have been important habitat patches in the dry periods, maintaining connectivity among populations, facilitating colonization of new habitats, and promoting the expansion of populations from central to southeastern Brazil.
One caveat to consider is that our divergence time estimates are based on calibrated mutation rates for two mitochondrial genes derived from other frog species. Our estimates under different evolutionary models were fairly consistent, but significantly more sensitive to assumptions of individual gene mutation rates (Table 7). In addition, the divergence time estimates from IM for node C are very different than the Tmrca estimates from BEAST. In part, this discordance is explained by the difference between divergence times (the time at which populations diverged) and Tmrca (gene divergence times); the first determines when populations were established, and the second determines when those populations become genetically structured (Wakeley & Hey 1997;Edwards & Beerli 2000). Populations may continue to experience gene flow following the founding of populations, or may have been isolated, but become panmictic again in Pleistocene refugia causing the two temporal estimates to diverge. Therefore, in the absence of fossil calibration, divergence times and Tmrca estimates should be considered an indication of relative divergence times, but they do not represent absolute dates.

Molecular diversity and historical demographic processes
As expected for species living in habitats that have historically undergone contractions and expansions, our mismatch analysis and ⁄ or Tajima's D show that both SE clade and CC populations expanded. We infer that this was likely a result of habitat expansion during that time. Likewise, IM results show that the ancestral CC ⁄ SE population was significantly smaller than the population sizes of those two daughter populations, a pattern consistent with demographic and geographical expansion from the ancestral range. In contrast, CG populations show the genetic signature of historical stability in population size. Therefore, older and more recent divergences in this species group resulted from very different evolutionary processes.
The patterns of diversification we uncovered indicate that landscape has a profound influence on the genetic structure of H. albopunctatus populations. The genetic divergence of the CG clade correlates qualitatively with the formation of the Chapada dos Guimarães. Likewise, populations in the SE clade are highly structured, with at least three well-supported monophyletic subclades. The southeast harbours a complex of mountain ranges parallel to the coast, with elevations reaching up to 2000 m that originated in the Tertiary (Safford 1999;Saadi et al. 2005) and are within the pre-settlement distribution of the Atlantic forest domain (Fig. 1). Although H. albopunctatus do not occur above 1000 m, the genetic structure among SE populations coincides with these highlands (subclade SE1 on the central plateau of São Paulo and Serra do Cipó mountain range in Minas Gerias, SE2 in the Serra de Paranapiacaba mountain range in southern São Paulo, extending into the southern mountain range, and SE3 restricted to the Serra da Mantiqueira range in northeastern São Paulo and southern Minas Gerais; Fig. 2). The break in southern São Paulo state between subclade SE2 and the other SE clades to the north is similar to phylogeographical breaks reported for other Atlantic forest taxa (e.g. Grazziotin et al. 2006;Brunes et al. 2010;Thomé et al. 2010). The mid-elevation habitat preferences of this species may limit gene flow among populations, promoting localized genetic structure in areas with complex topography. In contrast, in the Cerrado of central Brazil, topography is less accentuated and habitat availability may be more continuous, perhaps facilitating connectivity and gene flow among populations (Fig. 1). This hypothesis is corroborated by the fact that CC haplotypes show similar levels of haplotype and nucleotide diversity, yet they are consistently less geographically structured than populations in other regions. Only two well-supported haplotype groups are evident among the CC haplotypes: CC1, including haplotypes from populations in the central Brazilian plateau, and CC2, composed primarily of samples from the Serra da Mantiqueira, plus one haplotype from central São Paulo and one from Serra do Cipó , Minas Gerais state (Fig. 2). These admixed regions either represent contact zones between CC and SE lineages with introgression in particular populations, or they represent lineage sorting in populations or groups of populations that have not yet reached reciprocal monophyly.
Overall, H. albopunctatus populations exhibited high levels of mtDNA haplotypic diversity and varied levels of nucleotide diversity (Table 3). The SE and CC haplotype groups have high haplotype diversity and lower nucleotide diversity (Table 3), a signature of population bottlenecks followed by rapid population expansion and accumulation of mutations (Grant & Bowen 1998). The mixture of SE and CC haplotypes in some populations could represent secondary contact between differentiated lineages; however, the IM estimates of migration rates among populations were very low. Therefore, if secondary contact has occurred, it has not resulted in widespread admixture of these differentiated regions. The pattern of haplotype sharing is asymmetrical: individuals collected in the southeastern region occupied by the Atlantic forest domain more often have haplotypes with affinities to the CC, but the reverse is not common. This pattern is consistent with the evolution of the SE clade from CC source populations and is confirmed by the significant asymmetry in migration between CC and SE revealed by the IM analysis. Given the recency of this divergence, and the potential for gene flow between these adjacent domains, it is not surprising that the two regions share haplotypes.
Barriers to gene flow have clearly prompted regional differentiation at the largest scale, especially between the Chapada dos Guimarães clade and other populations. A pattern of high haplotypic and nucleotide diversity, as seen in the CG clade, is typical of species with relatively large populations and long evolutionary histories (Grant & Bowen 1998;Werneck et al. 2011). The individuals sampled from the single CG population are on average more genetically distinct among themselves than they are when compared with haplotypes sampled across the entire range of the species. The high diversity found in this population could be the result of mutation accumulation, given the long time this lineage has diversified in isolation (Fig. 2, Table 7). A similar scenario was suggested to explain high levels of mitochondrial divergence between stoats (Mustela erminea) from Britain and Ireland (Martínková et al. 2007). Further work in the western Cerrado is needed to identify geographical boundaries of this clade and the processes maintaining this high inter-individual variation.
In contrast to patterns in the mtDNA fragments, haplotypic and nucleotide diversity of nuclear markers were low. Even individuals from the Chapada dos Guimarães exhibited low nuclear diversity, generally sharing common nuclear alleles found both in the CC and SE regions (Fig. 4). Discrepancies between mitochondrial and nuclear markers are not uncommon in recently diverged populations or closely related species (e.g., Monsen & Blouin 2003;Brunes et al. 2010). Although almost no information is available on the substitution rates of amphibian nuclear genes (Stö ck et al. 2011), it is evident that their evolutionary rates are much lower compared with mtDNA (Hoegg et al. 2004). Also, demographic changes, such as bottlenecks or founder effects, will have a larger effect on mtDNA because of smaller effective population sizes and higher mutation rates (e. g., Fay & Wu 1999;Prugnolle & Meeus 2002). Finally, discordant patterns in mtDNA and nuclear genes may also result from male-biased dispersal, which often occurs in polygynous species (Vos et al. 2000;Prugnolle & Meeus 2002;Lampert et al. 2003). Hypsiboas albopunctatus is a highly polygynous species and males may be more prone to move from one site to another, searching for females (C. P. A. Prado, personal observations). However, sex-biased dispersal has not been measured in most amphibians and discrepancies between bi-parental and uni-parental markers may also be influenced by differences in sex ratio and fluctuating effective population size (Monsen & Blouin 2003). Thus, sex-biased dispersal remains a potential mechanism, but not the only explanation for the discordance we observed between genetic markers.

The role of Tertiary and Pleistocene events in the diversification of South American open habitats
High species diversity in the Neotropics has been attributed to three major biogeographical events: the uplift of South American mountains and plateaus and the emergence of the Isthmus of Panamá during the Tertiary (e.g. Vanzolini & Heyer 1985;Silva 1995;Nores 2004), and the more recent Pleistocene climatic fluctuations in tropical South America (Vanzolini & Williams 1981;Haffer 1997). A recent review of the timing of divergences in Neotropical taxa revealed that this dichotomy between older geographical events and more recent ones is an oversimplification: Neotropical species have originated continually since the late Eocene and the present day, and the biota is a mixture of species of different ages that diversified through various mechanisms (Rull 2008). In addition to the two tectonic events mentioned earlier, other potentially important Tertiary events include the uplift of the central Brazilian plateau and the Brazilian coastal mountain ranges with the formation of the adjacent depressions, alterations in river basins, marine introgressions, the origin of the savannalike habitat, and changes in forest distribution (Webb 1995;Saadi et al. 2005;Garda & Cannatella 2007). The deepest divergences in H. albopunctatus are relatively old (5.8 and 2.5 Ma) and thus may have been the result of these landscape changes during the Tertiary. Thus, our results corroborate similar findings of Tertiary events leading to vicariant diversification of open-area vertebrates in South America (Silva 1997;Porzecanski & Cracraft 2005;Werneck & Colli 2006;Maciel et al. 2010).
The nature and significance of Pleistocene climate cycles are still widely discussed in the literature (Colinvaux et al. 2001;Carnaval et al. 2009;Thomé et al. 2010). The Pleistocene Refugia Model (PRM) (Vanzolini & Williams 1981) is based on the premise that contraction and expansion of forests caused by climatic changes during the glacial and interglacial periods of the Pleistocene resulted in vicariant speciation. Recent molecular studies found speciation events concordant with the PRM (e.g. Ribas & Miyaki 2004;Almeida et al. 2007), whereas many others conclude that the Pleistocene climate changes may have contributed to population-level diversification, but not speciation (e.g. Eberhard & Bermingham 2004;Wü ster et al. 2005;Carnaval et al. 2009). Each species' response to changes in climate will depend in part on species-specific physiological and ecological requirements; thus, we would not expect to find a pattern that could explain diversification in the Neotropics as a whole. However, a pattern for the diversification of open habitats in South America is emerging from published phylogeographical studies of plants and animals (Ramos et al. 2007;Moraes et al. 2009). Molecular studies of typical Cerrado trees indicate the influence of Pleistocene climatic changes in population diversification, with lineages now present in the South and Southeast Brazil originating from eastern, central and western Cerrado (Collevatti et al. 2003;Ramos et al. 2007), a pattern similar to that we found for H. albopunctatus populations. A phylogeographical study of the cactophilic species Drosophila gouveai, which is associated with xeric vegetation from the Caatinga and enclaves in the Cerrado of Southeast Brazil, revealed signals of population expansion and diversification in the late to mid-Pleistocene (Moraes et al. 2009). Our phylogeographical results show similar genetic signature of Pleistocene habitat changes and confirm that these cycles can be just as significant in determining population structure in open areas as in forest-adapted species.
In summary, the treefrog H. albopunctatus shows a pattern of surprisingly deep divergences, cryptic diversity, and regional population structure despite a widespread distribution and generalized ecology. The history of this species has been characterized by vicariant isolation of lineages, followed by population expansion into new regions, and genetic differentiation associated with migration into new expanded habitats. Our findings confirm that the Cerrado is not composed of continuous and barrier-free habitat and that the species living there have had a more complex evolutionary history than previously predicted (Werneck 2011). Historical changes in the distribution and amount of open habitats in Central Brazil contributed to overall levels of diversification in a number of savanna species , and diversity in the Cerrado may be higher than currently recognized (Valdujo 2011). Further fine-scale studies of savanna-adapted taxa will confirm the generality of the spatial divergence patterns and timing we found here and further elucidate the role of historical habitat changes in the diversification of South American biota. This collaboration stems from our combined interest in processes generating diversity in Brazilian anurans, and was part of C.P.'s postdoctoral research on phylogeography of Cerrado frogs. C.P. is interested in ecology and evolution of open-habitat Brazilian frogs. C.H.'s laboratory studies taxonomy, ecology, natural history and conservation of Neotropical frogs. K.Z.'s laboratory studies patterns of diversification and population genetics of reptiles and amphibians in the New World.

Supporting information
Additional supporting information may be found in the online version of this article. Please note: Wiley-Blackwell are not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.