INTRODUCTION

Rivers are one of the most prominent landscape features, acting as key promoters of diversification among aquatic (Oliveira et al. 2019), semiaquatic (Roberto et al. 2020), flying (Hayes and Sewlal 2004), and terrestrial organisms in the Neotropics (Coelho et al. 2022). A er the first realization that primate ranges in the Amazon are delimited by major rivers (Wallace 1852), many studies have evaluated the Riverine Barrier Hypothesis (RBH) using distributional and genetic data of terrestrial, semiaquatic, and flying organisms in the Neotropics and elsewhere (Kopuchian et al. 2020, Coelho et al. 2022, Janiak et al. 2022). Rivers also represent the boundaries of areas of endemism and species turnover in Amazonia and the Atlantic Forest (Silva et al. 2005, Carnaval et al. 2014), while isolating parts of the Cerrado and compartmentalizing the landscape into highlands and lowlands (Brasil and Alvarenga 1989). Along with the watercourse features (e.g. size, volume, and flow direction), the ecological and natural history idiosyncrasies of the taxa in question also help define effects of rivers on population structure.

Accordingly, the rise of intraspecific phylogeography saw the RBH as one of the first and most directly testable mechanisms to be examined using mitochondrial DNA (mtDNA) sequences (Gascon et al. 1998, 2000, Lougheed et al. 1999). Indeed, rivers have been shown to drive the differentiation of aquatic or semiaquatic fauna in taxa such as caimans (Muniz et al. 2017, Roberto et al. 2020), snakes (Brandley et al. 2010, Ukuwela et al. 2013), fish (Hubert et al. 2007, Santos et al. 2009, Dubut et al. 2012), insects (Price et al. 2010, Metcalfe et al. 2020), anurans (Fonseca et al. 2021), and birds (Choueri et al. 2017(Choueri et al. , om et al. 2020)). Although rivers have been frequently shown to act as hard and so barriers to gene flow (Pyron and Burbrink 2010), they have also presented no effect on the distribution of genetic diversity in other cases (Fluck et al. 2020). e degree of genetic structure on opposing river banks has been shown to vary throughout the course of rivers, with middle and lower portions of the huge Amazonian rivers, for example, impeding migration more across margins than upstream regions (Reis et al. 2019).

Besides acting as barriers for species and populations across opposing margins, river basins have been shown to facilitate gene flow in semiaquatic and water-dependent species (Lawson 2013, Fonseca et al. 2021). At the landscape level, hydrological drainages were shown to function as corridors of connectivity among populations within basins (Lawson 2013, Choueri et al. 2017), with migration across populations preferably occurring downstream within each river basin ( om et al. 2020( om et al. , Fonseca et al. 2021)). is alternative role of rivers was recently named the Riverine ruway Hypothesis (RTH), which posits that rivers should act as facilitators of gene flow in species associated with aquatic environments (Fonseca et al. 2021), which include fish (Hubert et al. 2007) and even birds ( om et al. 2020).

Still, most evolutionary studies on anurans in the Neotropical region have focused on rivers as vicariant barriers (Kaefer et al. 2012, Moraes et al. 2016, Maia et al. 2017, Godinho and da Silva 2018). Furthermore, few studies have examined aquatic or semiaquatic anurans, even though some of these species are strongly associated with rivers and floodplains (e.g. Pseudae, Lithobates spp., Pipa spp. : Trueb and Cannatella 1986, Hillis and Wilcox 2005, Garda et al. 2010). us, such semiaquatic species with strong ecological bounds associated with floodplains are expected to follow river histories and be restricted to specific basins.

e South American paradoxical frogs (Hylidae: Pseudinae) are morphologically adapted to live in river floodplains (Aguiar et al. 2007, Garda and Cannatella 2007, Garda et al. 2010). eir large larvae (up to 270 mm in total length, Santana et al. 2016), shrink up to five times their length during metamorphosis into an ordinary-sized adult frog (32.3-52.9 mm snout-vent length, Garda et al. 2010). Such habits and geographical distribution make these frogs ideal organisms to test hypotheses of several evolutionary pa erns, including biogeographical hypotheses regarding speciation of semiaquatic fauna. Gallardo (1961), for example, suggested a 'one basin-one species' hypothesis for the group, a pa ern usually found in strictly aquatic organisms. Recently, a landscape genetics study with Pseudis tocantins Caramaschi and Cruz, 1998, which is distributed throughout the floodplains of long and wide rivers in central Brazil, showed that gene flow occurs primarily within basins (Fonseca et al. 2021). Among the seven species of the genus, Pseudis bolbodactyla Lutz, 1925 occurs in the floodplains of the middle São Francisco River, in the upper Tocantins River, which is part of the upper Araguaia-Tocantins basin , and along the Paranaíba River, which is part of the upper Paraná basin (Garda et al. 2010), all of which are watersheds restricted to the Brazilian Shield.

Herein, we test the RTH (i.e. rivers as facilitators to gene flow) and the 'one basin-one species' hypothesis using P. bolbodactyla as an organism model. If both hypotheses are corroborated, we expect to see each basin harbouring a distinct lineage with li le to no gene flow among basins. We also explore phylogeographical pa erns and the underlying processes that drove the current levels of genetic variation within P. bolbodactyla. To achieve these goals, we used phylogeographic and phylogenetic analyses to perform population assignment tests, estimate evolutionary relationships and historical demographies, and reconstruct evolutionary history using a Bayesian-based method and testing competing evolutionary models using approximate Bayesian computation (ABC).

MATERIAL AND METHODS

Sample collection and sequencing

We collected 166 tissue samples of P. bolbodactyla from 13 localities (Fig. 1; Supporting Information, Table S1; collecting permits #132/2005 CGFAU/LIC and SISBIO 26157). We used the sister species Pseudis paradoxa (KP149190) as the outgroup to P. bolbodactyla (Garda and Cannatella 2007). Voucher specimens from our field trips are housed in Coleção Herpetológica da Universidade de Brasília, Brasília, Brazil (CHUNB; Appendix 1). Whole genomic DNA was taken from muscle or liver tissues using the Qiagen DNeasy Blood and Tissue Kit (Qiagen). We followed the standard protocol of polymerase chain reaction (PCR) for DNA amplification for one mitochondrial marker [cytochrome c oxidase subunit 1 (COI)] and two nuclear markers [recombination activating gene 1 ( G-1) and proopiomelanocortin (POMC)]. Information on primers used in this study is presented in Table S2. PCR products were cleaned using Exo-Sap-IT (USB Corp.) and sequenced in both directions on a Beckman-Coulter CEQ-8000 automated sequencer. We performed all sequence alignments in Geneious ® v.9.1.8 using the MUSCLE alignment algorithm. We used PHASE v.2.1.1 (Stephens and Donnelly 2003) to determine the most probable pair of alleles for both nuclear loci. We ran PHASE with default parameters for 100 iterations, a thinning interval of 1, and a burn-in of 100. To check for consistency between runs, we repeated each PHASE analysis five times. e most appropriate model of nucleotide substitution for the alignments was determined with jModeltest (Darriba et al. 2012) using the Bayesian Information Criterion (BIC).

Population assignment and genetic diversity

We performed the Bayesian clustering method implemented in the R package Geneland v.4.0.4 (Guillot et al. 2005) to detect geographical discontinuities in both mitochondrial and nuclear loci. is method considers the spatial coordinates of each individual sampled and distributes them into K clusters by minimizing deviations from Hardy-Weinberg equilibrium and gametic phase disequilibrium within the groups (Guillot et al. 2005). We used the uncorrelated allele frequency model and evaluated support for 1-5 populations with 100 iterations and a burn-in of 1000. Every 100th observation was sampled to reduce autocorrelation. To estimate K, we simulated a fixed value of K using the above parameters to determine population membership and generate posterior probability maps.

Genetic structure of Pseudis bolbodactyla • 3

Based on the lineages delimited according to the method above, haplotype diversity (Hd) and nucleotide diversity (π) were calculated for each group using DNAsp v.5 (Librado and Rozas 2009). We also performed an analysis of molecular variance (AMOVA, Excoffier et al. 1992) partitioning our samples into hierarchical components among lineages to assess how much genetic variation was present within and between lineages. Population substructure was quantified with the fixation index F ST in the program Arlequin (Librado and Rozas 2009) a er 10 000 iterations. Finally, we explored the relationships among haplotypes of each locus in POPART (Leigh et al. 2015) through the median-joining network method. We identified each lineage using different colours in the haplotype network. We also calculated sequence divergence (uncorrected p-distance) among lineages using MEGA X (Kumar et al. 2018) to assess the degree of genetic divergence.

We estimated a Bayesian gene tree for the mtDNA loci (COI) using BEAST v.2.6 (Bouckaert et al. 2019). We used a Yule speciation prior, implemented a strict clock rate of 0.00957 per Myr (Crawford 2003), and ran the analysis for 50 million generations sampling every 1000 generations. We used Tracer v.1.7.1 (Rambaut et al. 2018) to assess effective sample sizes (ESS) of estimated parameters and stationarity, ensuring that all ESS values of all parameters were above 200 (Rambaut et al. 2018).

Dated molecular phylogeny

e three lineages found in Geneland (herea er referred to as AT, SF, and PR; see Results below) were used as terminal taxa for the species tree estimation in *BEAST (Heled and Drummond 2010) using COI, POMC, and G-1 unliked with no partitions, and an evolutionary rate for COI in BEAST v.2.6.3. To infer the timing of lineages divergence we used the mtDNA COI mutation rate of 0.957% per lineage per Myr (Crawford 2003). We calibrated the species tree using this mutation rate because there are no fossils available for calibration. We ran the analysis with 200 million generations, sampling every 10 000 generations. Stationarity was determined by visually inspecting trace plots and ensuring that all ESS values were above 200 in Tracer v.1.7.1. e first 10% of sampled genealogies were discarded as burn-in, and the most credible clade was inferred with TreeAnnotator v.2.6.3 (Bouckaert et al. 2019). We used DENSITREE (Bouckaert 2010) to produce a tree cloud from sampled trees.

Historical demography and diversification scenarios

We evaluated past changes in the effective population size (Ne) of each lineage using the Bayesian Skyline Plot (BSP) method implemented in BEAST v.2.6.3. For the BSP we used an average substitution rate for vertebrate mtDNA of 0.957% per lineage per Myr (Crawford 2003) and the same substitution model for 4 • Santana et al.

COI. Each analysis of BSP was evaluated with a linear substitution model, five groups, 50 million generations, and 10% burn-in. We checked for stationarity by visually inspecting trace plots and ensuring that all values for ESS were above 200 in Tracer v.1.7.1. We took the values from Ne over time from Tracer v.1.7.1 and built the BSP curve in Microso Excel.

We used ABC to simulate and calculate the posterior probability of competing diversification scenarios. We specified a total of five scenarios with stable population sizes (Fig. 2): (i) no migration among lineages; (ii) gene flow from PR to the ancestor of lineages AT and SF; (iii) gene flow between lineages AT and SF; (iv) gene flow from PR to the ancestor of lineages AT and SF, and gene flow from SF to AT; and (v) gene flow from PR to the ancestor of lineages AT and SF, and gene flow between lineages AT and SF. We constructed these models based on the results of previous analyses and on the biogeographical history of the Paraná Basin. We built diversification scenarios that mirrored our empirical datasets given the number of genetic markers, marker inheritance, number of individuals per marker and lineage, and sequence length. All priors were sampled from a uniform distribution and the value of each prior is available in Supporting Information Table S3. We used the results recovered by *BEAST analysis to set the evolutionary relationships among lineages and both divergence time priors (T 1 and T 2 ). en, we used the R package PipeMaster (Gehara et al. 2020) for genetic simulations. We performed 100 000 simulations under each demographic scenario and calculated a total of 26 genetic summary statistics from each individual simulation: number of segregating sites (S), nucleotide diversity (π), haplotypic diversity, Tajima's D, Fu and Li's D and F statistics, and fixation index (F ST ). Summary statistics were calculated combining all lineages and for each lineage.

e only exception was F ST , which was calculated considering only individual lineages. Next, we used the function postpr implemented in the R package abc (Csilléry et al. 2010) to calculate the posterior probability of each diversification scenario given the empirical dataset. To compare our demographic scenarios, we used rejection and mnlogistic methods implemented in the abc R package with tolerances of 0.001 and 0.01, respectively. We evaluated the performance of our model by generating a confusion matrix using the cv4abc function in the abc R package using the mnlogistic method with a tolerance of 0.01. We also used a principal components analysis (PCA) to check if simulated datasets produced summary statistics similar to our observed dataset. Finally, we estimated the posterior probability of divergence times and migration rates under the best model.

RESULTS

Population assignment, genetic diversity and species tree

We successfully amplified gene fragments from 154 samples.

e mtDNA gene tree topology recovered three lineages corresponding to the three hydrographic basins: Paraná (PR), Araguaia-Tocantins (AT), and São Francisco (SF) River basins. Overall, these lineages were inferred with strong support [posterior probability (pp) > 0.99; Fig. 3]. However, two individuals from the São Francisco basin (CHUNB42877 and CHUNB42885) were more closely related to individuals from the Araguaia-Tocantins basin. Using all three loci, the Geneland analysis detected the same three lineages (K = 3; Fig. 4), revealing a clear geographical population structure matching hydrographic basins. Genetic distances among populations ranged from 3.3% to 7.7% (Table 1). e mtDNA COI haplotype network revealed high haplotype diversity, with no shared haplotypes among lineages (Fig. 5). e SF and PR lineages present one dominant haplotype each in the POMC haplotype network, while the PR lineage is more diverse. e G-1 haplotype network has one central and most frequent haplotype, which is shared among lineages. According to AMOVA (Table 2), most genetic variation in COI occurs among lineages (F ST = 0.82; P < 0.00001). e nuclear genes presented contrasting results. We found a high index of fixation for POMC (F ST = 0.76033; P < 0.00001) and a low index for G-1 (F ST = 0.03394; P = 0.03089).

e species tree reconstructed in *BEAST showed high support for AT, SF, and PR lineages (P = 1) and estimated the divergence time between the PR lineage and SF-AT around 3.08 Mya [95% highest posterior density (HPD) = 1.56-4.63], during the Plio-Pleistocene (Fig. 6). *BEAST estimated a more recent separation between SF and AT dated to 1.21 Mya (95% HPD = 0.53-1.95), during the Pleistocene.

Historical demography, migration, and diversification scenarios

e BSP for COI showed no evidence of demographic fluctuations for the SF and AT lineages, while PR showed a moderate expansion over the past 200,000 years, albeit within the 95% confidence interval (10× expansion around 0.4-0.2 Mya, Fig. 7). e AT and SF lineages have small effective sample sizes (Ne) in comparison with PR, which is almost 10 times larger.

e ABC analysis recovered model 5 (Fig. 2; Table 3) as the best diversification scenario, with a probability of 0.86 and 0.99 for rejection and mnlogistic methods, respectively (Table 3).

is model predicts a past migration from PR to the ancestor of AT-SF, and current bidirectional gene flow between AT and SF (Fig. 2). A rate of one migrant per generation was inferred for both gene flow from PR to the ancestor of AT-SF (95% HPD: 0.02-1.98) and for the bidirectional migration between AT and SF (95% HPD: 0.09-1.96). ABC analysis had an overall accuracy of 72%, with the best model reaching an accuracy of 92% (Supporting Information, Fig. S1). Our simulated datasets produced summary statistics that were in the range observed in our dataset (Fig. S2).

DISCUSSION

Pseudis bolbodactyla is composed of three deeply diverging lineages (AT, SF, and PR) associated with three unique hydrographic basins. We detected li le current gene flow between the AT and SF lineages and high genetic differentiation among hydrographic basins (except for G-1). is matches the predictions of the RTH and the 'one basin-one species' hypothesis, in which populations or species associated with aquatic environments are genetically structured by hydrographic basins. While rivers facilitate within-basin migration (RTH), populations become isolated in each basin (one basin-one species hypothesis). Both mechanisms therefore simultaneously promote diversification among but intensify gene flow within basins, hence accounting for the pa ern observed. is pa ern was also found for P. tocantins in the Tocantins-Araguaia River system (Fonseca et al. 2021). ese two rivers run in parallel northward for hundreds of kilometres before merging, and although separated only by moderately high mountains (up to 500 m), gene flow occurs within the basin and primarily downstream (Fonseca et al. 2021). While seemingly intuitive, the role of rivers in structuring populations and favouring gene flow is rarely tested explicitly for semiaquatic and terrestrial animals, with the more common finding showing rivers as barriers (Fonseca et al. 2021).

The lineages of P. bolbodactyla diverged during the Plio-Pleistocene (PR and AT-SF) and Pleistocene (AT and SF), periods marked by orogenic events in central Brazil that led to the present configuration of hydrographic basins (Lundberg et al. 1998). Pleistocene climatic fluctuations may have also influenced divergence, but in different ways for each basin. The Araguaia-Tocantins basin attained is current configuration during the end of the Pliocene and throughout the Pleistocene (Brasil andAlvarenga 1989, Del'Arco andBezerra 1989), allowing the divergence of AT and SF from the shared ancestor with PR under unidirectional migration from PR. Confinement to and migration among hydrographic basins in central Brazil were probably influenced by the upli of the central Brazilian plateau and related events of headwater capture (Aquino and Colli 2017). For example, AT has a lower baseline compared to the PR and SF basins, making it more prone to capture headwaters of neighbouring river basins (Saadi 1993, Aquino andColli 2017). Indeed, the fish fauna in PR is older and more structured, probably having colonized the upper SF and AT basins in headwater capture events (Aquino and Colli 2017). Our PR lineage is also older, and we detected ancient unidirectional migration from PR to the ancestor of AT-SF (model 5; Fig. 2), which suggests that headwater capture events must have played similar roles in the current pa ern of genetic diversity and distribution in P. bolbodactyla. Our AT and SF lineages diverged more recently with low bidirectional migration until the present day (one migrant per generation; see ABC results) but remained limited to the upper Tocantins River and São Francisco Rivers, respectively. is restricted distribution of AT and SF populations has two possible explanations. First, the middle and lower parts of the São Francisco River run across the Brazilian semiarid Caatinga region, which could limit dispersion into such areas probably due to prolonged and severe droughts reducing the occurrence and duration of suitable water bodies. Second, the Paranã Valley in the upper Tocantins River is isolated from downstream sites by a narrow canyon, which slows the river and forms a large wetland in northeastern Goiás State.

ere is limited information on the fauna shared by these three hydrographic basins other than fishes (Aquino and Colli 2017).

e Cerrado geomorphological compartmentalization resulting from early Miocene tectonic activities has been implicated in lineage diversification and speciation in this region (Guarnizo et al. 2016). is activity led to upli of the Brazilian Shield, which subsequently started a long process of erosion that shaped the landscape into ancient plateaus, dominated by savanna-like vegetation, and younger valleys with more heterogeneous forest assemblages (Colli 2005). eoretically, this dynamic history could have resulted in distinct lineages forming in plateaus and valleys (Werneck et al. 2011). However, different diversification pa erns have been inferred from phylogeographical assessments of the Cerrado biota, with clades changing geographically in two directionality pa erns: northwest-southwest (for Cerrado endemics) and southwest-northeast (in groups widespread along Genetic structure of Pseudis bolbodactyla • 9 the diagonal of open formations, encompassing the Caatinga, Cerrado, and Chaco biomes, Guarnizo et al. 2016). Our results show a different pa ern for a Cerrado endemic species, with a southeast-northeast directionality shaped by hydrographic basins and demographically influenced by climatic fluctuations during the Pleistocene. Plants and animals in the diagonal of open formations show idiosyncratic responses to Pleistocene climatic cycles, but most have signs of synchronous population changes (Gehara et al. 2017, Bonatelli et al. 2022). ere is no evidence of demographic fluctuation for the SF and AT populations, while PR shows a moderate expansion over the past 200 000 years. e Paraná basin runs southward, and some areas currently occupied by the species probably experienced colder climates and/or stronger influences of tropical forests during the Pleistocene. As climate changed and forests receded, adequate habitats became available for these tropical, floodplain-dwelling frog species. e timing of this putative expansion for the PR population matches those observed for other open-area species, such as Caatinga dry forest frogs and lizards (Gehara et al. 2017).

Such a combined effect of demographic histories in response to climatic fluctuations and the Brazilian plateau upli is common among many vertebrates in the region (Maciel et al. 2010, Prado et al. 2012, Oliveira et al. 2018). In cane toads [Rhinella marina (Linneus 1758) species group], for example, two events of diversification were caused by upli of the Brazilian Shield during the Miocene, which split them into two main lineages (north and south of the Brazilian Shield, Maciel et al. 2010). e Cerrado treefrog Boana albopunctata (Spix 1824) shows the same diversification process, with central and southeast Cerrado clades (Prado et al. 2012). Within the fossorial but still widely distributed frog Dermatonotus muelleri (Boe ger, 1885), the Central Brazilian plateau acted as a barrier between two distinct lineages across the open formations of South America (Oliveira et al. 2018). Pa erns recovered for these species, as well as the timing of events, are similar to what we recovered herein for P. bolbodactyla, corroborating Miocene orogenic events as drivers of anuran diversification in Central Brazil. e recent diversification in the Pleistocene for AT and SF probably reflects a concordant geomorphological history between these two drainages, which diverged later than PR drainages (Brasil and Alvarenga 1989, Del' Arco and Bezerra 1989, Turche o-Zolet et al. 2013).

e ancient separation of the PR basin (Gray et al. 1985, Po er 1997) and its distinct geomorphological history (Gray et al. 1985, Campos and Dardenne 1997, Lundberg et al. 1998) likely influenced the formation of this lineage.

Taxonomic implications of the one basin-one species hypothesis

Recent taxonomic revisions of paradoxical frogs have compared all species and many populations using morphological, morphometric, larval and bioacoustic characters, and molecular data (Garda et al. 2010, Santana et al. 2013, 2016). Despite populations of some species being geographically well structured, many cannot be currently distinguished without molecular data. For example, Lysapsus bolivianus Gallardo 1961 and L. limellum Cope 1862 are virtually identical based on larval, morphological, and bioacoustic characters (Santana et al. 2013). Likewise, the sister species P. paradoxa (Linnaeus 1858) and P. platensis Gallardo 1961 and populations of P. bolbodactyla are morphologically and bioacoustically cryptic (Santana et al. 2016).

In this study, we found that P. bolbodactyla is a complex formed by three distinct lineages, each restricted to a distinct hydrographic basin in central Brazil. Such results corroborate the predictions of the RTH and indicate that its scenarios must be considered as well when studying the landscape genetics of Neotropical species, especially semiaquatic ones. ese lineages might correspond to undescribed species, a hypothesis that can be further tested using denser genetic and morphological sampling. As two samples from SF were nested within AT, suggesting gene flow or incomplete lineage sorting, future studies should incorporate genomic data to investigate speciation with gene flow while testing for other factors such as isolation by distance and isolation by resistance (plateaus vs. floodplains). Nonetheless, we suggest that large genetic distances, distinct demographic histories, and the likelihood of allopatric speciation may be similar to those seen in other aquatic and semiaquatic tetrapod fauna such as frogs (Pipa, Lithobates, Lysapsus), cecilians (Typhlonectes, Atretochoana), lizards (Dracaena, Crocodilurus, Potamites, Neusticurus), and snakes (Helicops, Erythrolamprus).

SUPPLE MEN TARY DATA

Supplementary data are available at Biological Journal of the Linnean Society online.