Biogeography of curimatid fishes reveals multiple lowland–upland river transitions and differential diversification in the Neotropics (Teleostei, Curimatidae)

Abstract The Neotropics harbors a megadiverse ichthyofauna comprising over 6300 species with approximately 80% in just three taxonomic orders within the clade Characiphysi. This highly diverse group has evolved in tropical South America over tens to hundreds of millions of years influenced mostly by re‐arrangements of river drainages in lowland and upland systems. In this study, we investigate patterns of spatial diversification in Neotropical freshwater fishes in the family Curimatidae, a species‐rich clade of the order Characiformes. Specifically, we examined ancestral areas, dispersal events, and shifts in species richness using spatially explicit biogeographic and macroevolutionary models to determine whether lowlands–uplands serve as museums or cradles of diversification for curimatids. We used fossil information to estimate divergence times in BEAST, multiple time‐stratified models of geographic range evolution in BioGeoBEARS, and alternative models of geographic state‐dependent speciation and extinction in GeoHiSSE. Our results suggest that the most recent common ancestor of curimatids originated in the Late Cretaceous likely in lowland paleodrainages of northwestern South America. Dispersals from lowland to upland river basins of the Brazilian and Guiana shields occurred repeatedly across independently evolving lineages in the Cenozoic. Colonization of upland drainages was often coupled with increased rates of net diversification in species‐rich genera such as Cyphocharax and Steindachnerina. Our findings demonstrate that colonization of novel aquatic environments at higher elevations is associated with an increased rate of diversification, although this pattern is clade‐dependent and driven mostly by allopatric speciation. Curimatids reinforce an emerging perspective that Amazonian lowlands act as a museum by accumulating species along time, whereas the transitions to uplands stimulate higher net diversification rates and lineage diversification.

As with all continental biotas, variation in both climate and landscapes affects the formation of NFF assemblages of lowland and upland river basins. Lowlands are defined as low-elevation (<250-300 m) wetlands, floodplains, lakes, and large river channels in the Amazon, Orinoco, and Paraguay basins, including large tributaries such as Madeira, Negro, Purus, Putumayo, Solimões, and Ucayali rivers. Uplands are defined as higher elevation (>250-300 m) basins associated with the Andean cordilleras (e.g., Magdalena-Cauca, Upper Meta) or cratonic highlands covering most of central and southeastern Brazil (the Brazilian Shield; e.g., Paraná, São Francisco, Upper Tapajós, Araguaia-Tocantins, Upper Xingu rivers) and portions of Venezuela and the Guianas (the Guiana Shield; e.g., Branco, Corantijn, Essequibo, Marowijne, and Oyapock rivers) . These elevation zones exhibit distinct geomorphological and environmental profiles and contrasting biogeographic patterns. Compared with upland basins, lowland basins are more centrally located on the continental platform, have larger rivers and floodplains with a greater total volume of aquatic habitat, and flow across more erodible alluvial substrates with more rapidly changing river courses. Lowland basins of northern South America (e.g., Amazon and Orinoco basins) also exhibit much higher species densities (i.e., local alpha diversity) and lower percent species endemism (i.e., geographic beta diversity) than do upland basins of the continental periphery (Albert et al., , 2020. The evolutionary origins of the remarkable NFF diversity remain incompletely understood, with several hypotheses linking mechanisms of diversification to regional geomorphological and climate histories Albert et al., 2020;Lundberg et al., 1998). One hypothesis states that the older cratonic Brazilian and Guiana shields are more stable geologically and environmentally, and therefore represent an "area of origin" for the NFF fauna (Eigenmann & Allen, 1942;Albert & Carvalho, 2011: fig. 7.3). Another hypothesis states that rivers in the large lowland basins undergo faster rates of tributary exchanges due to river capture and sea-level fluctuations, and therefore experience greater rates of both speciation and extinction than upland basins Lima & Ribeiro, 2011). Thus, lowland basins have been hypothesized to generate species more rapidly, acting as evolutionary cradles of diversity, while upland basins have been hypothesized to isolate species and protect them from extinction acting more as evolutionary museums of diversity (see Lundberg et al., 1998).
An evolutionary museum is a region of net species accumulation, in which regional extinction rates are lower than the combined rates of in situ speciation and immigration from adjacent regions (Stebbins, 1974). Conversely, an evolutionary cradle is a region of net species overproduction, in which regional speciation rates exceed extinction rates, resulting in an increase of species richness and biogeographic range expansion (i.e., dispersal) to adjacent regions ( Figure 1) Rangel et al., 2018;Stebbins, 1974). Although these simple macroevolutionary models do not fully summarize the complexity of most biogeographic regions, the museum-cradle paradigm is widely used in contemporary analyses of tropical diversification (Azevedo et al., 2020;Cássia-Silva et al., 2020;Dagallier et al., 2020;Matos-Maraví et al., 2019;Melo et al., 2021;Meseguer et al., 2020;Rangel et al., 2018). Some caveats that apply to these models are that the taxa of a given region are not expected to be monophyletic with respect to that of adjacent regions ( Figure 1)  and that a region may simultaneously serve as a cradle for some taxa and a museum for others (McKenna & Farrell, 2006;Moreau & Bell, 2013). Although the original museum/cradle concept of Stebbins (1974) was qualitative and descriptive, these terms are now used as shorthand for quantitative models that compare rates of net diversification among regions (e.g., Jablonski et al., 2006;Moreau & Bell, 2013). Although these macroevolutionary models have been tested in other tropical groups (Gaston & Blackburn, 1996;McKenna & Farrell, 2006;Moreau & Bell, 2013), none have yet been tested explicitly in the NFF fauna (see , but see Melo et al., 2021. Curimatid fishes, the branquinhas, saguirus, or toothless characiforms (Teleostei: Ostariophysi), are well suited to studies documenting the spatial and temporal contexts of evolutionary diversification. The Curimatidae is a broadly distributed clade of riverine-adapted characiform fishes that inhabits both lowland and upland river basins of tropical South and lower Central America (Vari, 1988). Members of the clade are detritivorous and perform short lateral migrations (Fernandes, 1997), with some species attaining larger body [up to 32 cm standard length in Curimata mivartii] and population sizes, thereby constituting an important component of Amazonian riverine food webs and human fisheries (Araujo- Lima et al., 1986;Araujo-Lima & Ruffino, 2003;Lowe-McConnell, 1975). Curimatids are important parts of the fish faunas in the lowland rivers of the Amazon, Orinoco, and La Plata basins, upland rivers of the Guiana and Brazilian shields (e.g., Essequibo, Tocantins, São Francisco, and Paraná basins), and swiftly flowing rivers west of the Andes (e.g., Maracaibo, Magdalena, Atrato and Pacific coastal basins from northern Peru to Costa Rica) (Frable, 2017;Vari, 1988). The alpha taxonomy and higher-level systematics of curimatids have attracted substantial attention (e.g., Vari et al., 2012;Melo, 2017Melo, , 2020Melo & Oliveira, 2017;Bortolo et al., 2018) using both morphological and molecular evidence (Dillman et al., 2016;Dorini et al., 2020;Melo et al., 2016Melo et al., , 2018Vari, 1989). The family is represented by 117 valid species in eight extant genera, representing one of the ten most species-rich families in the Amazon basin (Dagosta & de Pinna, 2019;Fricke et al., 2021;Vari, 1989).
Based on a morphological phylogeny, Vari (1988) investigated the historical biogeography of curimatids and his findings supported the delimitation of ten major bioregions in South America and the hypothesis that post-Andean cladogenesis had moderate effect on the curimatid diversification. Vari (1988) further discussed the close relationships among fish assemblages of the Amazon, Orinoco, and Guianas, and the mixed biota for the São Francisco basin and rivers of northeastern Brazil. More recently, a molecular phylogeny refined the relationships of curimatids including ca. 67% of all species diversity , resulting in a robust evolutionary framework to reassess Vari's hypotheses and to investigate spatial macroevolutionary scenarios for curimatid fishes. Combining this phylogeny  with an extensive spatial database, we estimate temporal diversification and geographic range evolution of curimatids aiming to elucidate F I G U R E 1 Alternative macroevolutionary models of species inhabiting lowland and upland river basins of tropical South America (from ; topologies reproduced with permission): (a) vicariance into lowland and upland descendants from geographically widespread ancestors; (b) lowlands as cradles of diversity with high diversification rates and uplands as museums with low diversification rates; and (c) lowlands as museums and uplands as cradles. Note other models are possible in which rates of speciation and extinction are decoupled (Louca & Pennell, 2020) whether lowland/upland river basins act as museum or cradles of lineage diversification in the Neotropics. Our findings on curimatid fishes support an emerging perspective of Amazonian lowlands as both a cradle and a museum of diversification, accumulating species over tens of millions of years, and serving the primary source of Neotropical biodiversity (Antonelli, Ariza, et al., 2018;, while the shield uplands act more often as a cradle of species diversification.

| Divergence time estimates
This study used as the phylogenetic framework a published multilocus dataset containing three mitochondrial loci, the 16S rRNA, cytochrome C oxidase subunit I (COI), and cytochrome B (Cytb), and three nuclear coding genes, the myosin heavy chain 6 gene (Myh6), recombination activating gene 1 (Rag1), and recombination activating gene 2 (Rag2), with a total of 5358 characters and 151 terminals . The original matrix was reduced to 142 terminals by excluding nine outgroup taxa because our aim was to estimate divergence times only. Voucher information and GenBank accession numbers are available in Table S1, and the additional locality information of taxa is available in Melo et al. (2018). PartitionFinder v2.1.1 (Lanfear et al., 2012) was used to infer the optimal partitioning scheme and model of DNA evolution for partitions, in which the non-coding 16S rRNA was set as one partition and the remaining five coding genes were subdivided by codon position. This analysis used the greedy algorithm and the Akaike information criterion (AIC) assuming the models available in BEAST v1.8.2 (Drummond et al., 2012) (Supporting Information). The best maximum-likelihood (ML) tree was obtained using the GTRGAMMA model in RAxML v8 (Stamatakis, 2014). An uncorrelated relaxed molecular clock with lognormal distribution was applied in BEAST v1.8.2 (Drummond et al., 2012), and the best ML tree was fixed as the topology with BEAST estimating only node ages.
This species shares with congeners modifications in the laterosensory canals of the fourth and fifth infraorbitals (Malabarba, 1996;Vari, 1989). †Cyphocharax mosesi has been tentatively hypothesized as close related to the extant species C. gilbert, C. modestus, and/or C. santacatarinae (Malabarba, 1996). The most recent molecular phylogeny of Curimatidae suggests these three extant species belong to the C. gilbert clade, which further includes C. corumbae, C. naegelii, C. platanus, C. spilotus, and C. voga . Thus, we used stratigraphic information of †Cyphocharax mosesi to constrain age estimates of the most recent common ancestor (MRCA) of C. gilbert clade using a lognormal prior (offset = 23.8 Ma; mean = 5.0; stdev = 7.5).
The second calibration is a constraint on the root of the tree selected to match the timing of the split between Curimatidae and Chilodontidae, as estimated from the molecular analysis of Characiformes . This analysis included 356 taxa and six fossil calibrations, and the time-calibrated phylogeny revealed the split of Chilodontidae and Curimatidae at approximately 54.9 Ma (64.2-46.2 Ma, 95% highest posterior density; HPD) ; therefore, we assigned a normally distributed prior on the root of our phylogeny (offset = 54.9 Ma; stdev = 7.5). The other available curimatid fossil †Plesiocurimata alvarengai (Figueiredo & Costa-Carvalho, 1999) is controversial because of the incompleteness of informative characters preventing any phylogenetic placement of the fossil species within the family (Figueiredo & Costa-Carvalho, 1999). The BEAST analysis was partitioned under the models of molecular evolution as determined by PartitionFinder, and a birth-death model was applied for the tree topology with branch lengths. BEAST ran 600 million generations, sampling frequencies at every 60,000th generation, which  Vari (1988) proposed 10 regions of endemism for curimatids. We tested this hypothesis with a dataset comprising 5362 GPS coordinates from 103 curimatid species assembled from museum collections, publications, and metadata repositories such as SpeciesLink (http://www.splink.org.br). We used the R package "Biogeo" (Robertson et al., 2016) to remove duplicates and records with obvious georeferenced errors. This procedure excluded occurrences in the ocean or outside the Neotropical region, those without country names, coordinates with zero latitude or longitude, and coordinates annotated on coarse-scale grid without decimal precision ( Figure 2a).

| Biogeographic regions
We used this assembled spatial dataset to delimit bioregions using bipartite network clustering in the Infomap Bioregions (Edler et al., 2016). This program uses an adaptive resolution of grid cells when sampling effort is uneven to generate bipartite networks (i.e., species vs. grid cells) clustered using the Infomap clustering algorithm.

| Biogeographic estimates
We used the R package "BioGeoBEARS" (Matzke, 2013) to compare the likelihoods of six time-stratified biogeographic models. The time intervals represented six epochs, from Late Cretaceous to Holocene.
The epochs Pliocene, Pleistocene, and Holocene were grouped into one interval. Constraints on area connectivity (i.e., areas-allowed) were similar across all six epochs, except in the last time slice where range expansions were not allowed between the cis-and trans-Andean river basins. Each model incorporates a combination of three analytical procedures and the inclusion or not of the distance-based parameter rate "x" (Van Dam & Matzke, 2016). The dispersal matrix was built using estimates of centroid points and Euclidian distances between adjacent bioregions and then rescaled by the shortest measured distance. The analytical procedures encompassed three biogeographic models: DEC model (Ree & Smith, 2008) and the probabilistic implementations of the programs DIVA (Matzke, 2014;Ronquist, 1997) and BayArea (Landis et al., 2013;Matzke, 2014).
Each analytical procedure allows an exclusive cladogenetic process  Note: Models with ΔAIC scores less than 2-log likelihoods are equally supported as the best-fitting model. Each model comprised a unique combination of ML analytical procedures (DEC, DIVA, BayArea) and a distance-based dispersal parameter rate (x).
Abbreviation: AICc, Akaike information criterion with corrections for sample sizes; d, range expansion parameter; e, range contraction; LnL, log-likelihood; n, number of parameters; ΔAIC, delta AIC.

TA B L E 1
Model selection among six time-stratified biogeographic models unlikely to undergo a jump speciation process because their dispersal abilities are usually constrained by river connectivity between adjacent watersheds. None of the curimatid species included in the analyses inhabited more than five bioregions; therefore, a maximum of five combined bioregions was allowed at each node of the timecalibrated tree. Input files of the six time-stratified biogeographic models are available in the Supporting Information.

| Model selection
Estimates of AICc scores compared global likelihoods at the root node across the six time-stratified biogeographic models (Table 1).
Besides penalizing additional model parameters, AICc has corrections for small sample sizes. When sample sizes are large, AICc scores converge to standard AIC scores. The sample size was the number of terminal taxa in the time-calibrated tree. The six biogeographic models included either two or three parameters: range expansion rate (d), range contraction rate (e), and distance-based dispersal rate (x). The lowest AIC score indicated the best-fitting biogeographic model. ΔAIC scores less than 2-log likelihoods were considered equally supported (Burnham & Anderson, 2004). Two of the six time-stratified models were considered equally supported (DIVA and DIVA + x) (Table 1). Although the two models include DIVA as their analytical procedure, they differ by the distancebased dispersal rate parameter "x." Given that both models are statistically supported (Table 1), the simplest one with the fewer parameters (i.e., DIVA) was used to perform subsequent biogeographic analyses.

| Biogeographic Stochastic Mapping
Biogeographic stochastic mapping (BSM) is a simulation-based approach for estimating biogeographic scenarios conditioned on a timecalibrated tree, biogeographic model, and its parameters (Dupin et al., 2017). Specifically, it is applied to simulate multiple equiprobable scenarios of range evolution given the best-fitting model, therefore allowing statistical inferences of modes and counts of biogeographic events (Dupin et al., 2017). BSM scripts available in "BioGeoBEARS" (Matzke, 2013) were used to estimate modes and counts of biogeographic events in the curimatid tree. To calculate the probability of ancestral states for each node, we performed 100 BSM simulations under the simplest best-fitting model (i.e., DIVA; Table 1). Cladogenetic events were binned in two geographic categories: sympatric (within-basin) and vicariant (among-basin) events. These events were further categorized by altitudinal gradients and vicariant patterns.

| Diversification rates
We used the Geographic State-dependent Speciation and Extinction (GeoSSE/GeoHiSSE) models (Caetano et al., 2018) in the R packages "devtools," "hisse," and "diversitree" (Beaulieu & O'Meara, 2015;FitzJohn, 2012) to estimate net diversification rates, and to test the hypothesis of range-dependent diversification processes (i.e., lowland/upland). These models apply a three-state Markov Chain allowing for species categorization in distinct (e.g., A or B) or combined (e.g., A and B) geographic units. Here, we categorized each species in a matrix with: (1) geographic range only in lowlands, (2) geographic range only in uplands, or (3) geographic range in both lowlands and uplands. Model parameters were the following: (1) speciation rate in each geographic unit (sLowands or sUplands) and in both units (sLowland and sUplands), and (2) extinction rate in each geographic unit (xLowlands or xUplands), and dispersal rates between geographic units (dLowands = lowlands to uplands, dUplands = uplands to lowlands). The state-dependent models optimize orthogonal transformations of these variables such as the parameter of net turnover for the widespread range (e.g., 01) and extinction fractions of endemics (e.g., 0 or 1). There is no "extinction" parameter associated with widespread lineages (see details in Caetano et al., 2018).
The chronogram and geographic range matrix were used as the input data for testing four models-(M1) a range-independent diversification process without hidden states ( and (M4) range-dependent diversification process with hidden states and multiple rate states: This model assumes shifts in diversification across the tree, and these are dependent on the geographic ranges. The relative importance of each of these models to explain the variation observed in the data was evaluated by AIC weights.
Input data and scripts are available in the Supporting Information.

| Species richness per time
To and the literature for the unsampled taxa without information of geographic ranges (e.g., Vari, 1991Vari, , 1992Vari, , 2003. We estimated the number of species per million of year by dividing species counts by the mean crown age for seven major subclades of the Curimatidae: Psectrogaster, (6) Cyphocharax+Curimatella, and (7) Steindachnerina.  . We also excluded Cy. abramoides, Cy. nigripinnis, and Cy. multilineatus from this analysis due to their phylogenetic placement outside the Cyphocharax + Curimatella clade .

| Range evolution in river basins
The bipartite network clustering in the Infomap Bioregions delim-

| Sympatric and vicariant events
The number of within-basin and among-basin cladogenetic events is similar within curimatids ( Figure S1). This result may reflect the grouping of many lowland basins (e.g., Amazon basin) as a single bioregion rather than partitioned into multiple biogeographic units (see Dagosta & de Pinna, 2017). Within-basin speciation events were more common within lowland than upland basins, and among-basin speciation events are about equally distributed in lowland/lowland and lowland/upland river basins. Among within-basin events, the most common was, by far, within the Amazon lowlands, followed by Trans-Andean and Paraguay river basins ( Figure S1). The most common among-basin events were between Amazon/Paraguay followed by four transitions involving the Amazon basin: the vicariant events in Amazon/Araguaia-Tocantins, Amazon/Orinoco, Amazon/Trans-Andean, and Amazon/Guianas ( Figure S1).

| From lowlands to uplands
Using the best-fitting model (

| Diversification rates
The range-dependent diversification analysis (GeoSSE/GeoHiSSE) tested four distinct models: M1 (null model) assumes a homogeneous diversification rate across the tree independent of the ranges; M2 (GeoSSE) assumes a range-dependent diversification process without hidden states; M3 (GeoHiSSE) assumes shifts in diversification across the tree independent of the ranges; and M4 (GeoHiSSE) assumes a range-dependent diversification process with hidden and multiple rate states (Caetano et al., 2018). We computed the relative importance of each of the models to explain the variation observed in the data using AIC weights. Results indicate that the best-fitting model was M4 (Table 2), where variation in net diversification rates is correlated with geographic ranges (in our case, lowlands and uplands). The two best-fitting models (M2 and M4) account together for 0.999 of the AIC weight (0.243 and 0.756, respectively), assuming shifts in diversification rates, and vary in whether hidden states are allowed or not. The two other models treating diversification and range evolution as independent processes appeared with poorer fits and presented AIC weights below 0.0001 (Table 2). We then performed a marginal reconstruction for each of the models in the set, reconstructed the hidden states at the nodes of the phylogeny, and plotted the results in the topology ( Figure 5). Results indicate a range F I G U R E 3 Time-calibrated tree estimated in BEAST and ancestral range evolution of curimatids by bioregions estimated in BioGeoBEARS. The most likely ancestral state is colored in node circles, and less likely ancestral states are combined and shown in white color. Results show that the most recent common ancestor of all curimatids originated in the Late Cretaceous at around 78 Ma (90-66 Ma 95% HPD) in the lowlands of the Proto-Orinoco-Amazonas basin with multiple dispersal events to uplands of the Brazilian and Guiana shields of 0.0006 to 0.261 for net diversification rates with higher rates occurring mostly in clades with widespread taxa (i.e., taxa occurring in both lowland and upland basins), as denoted by the yellow branches with bright red outlines ( Figure 5). These are more evident in subclades of Cyphocharax and Steindachnerina with a few instances within Curimata and Pseudocurimata ( Figure 5).  (Table 3).

| Cretaceous origins and Cenozoic diversification
The time calibration and ancestral biogeographic estimations indicate an origin of curimatids in lowlands of the Proto-Orinoco-  (Malabarba, 1996;Melo et al., 2018).
The presence of †Cyphocharax mosesi, a fossil phylogenetically close to modern species near the Oligocene-Miocene boundary, suggests that at least some modern curimatid genera were already present by that time period. This perspective is also supported by the presence of †Plesiocurimata in the same Tremembé Formation, which may be phylogenetically closer to a clade of most derived taxa (Curimatella,
Most species-level diversification within curimatid genera is es-  (Table 3), although the actual diversity of this and other clades seems to be underestimated (Melo et al., 2016). We hypothesize that diversification of these clades may have been triggered by: (1) large river captures promoting dispersal from POA to adjacent upland drainages of the Guiana and Brazilian shields (Albert, Val, et al., 2018;Tagliacollo et al., 2015); (2) derived morphophysiological or behavioral adaptations that allowed members of colonizing populations access to novel food or habitat resources in upland basins (Silva et al., 2016), and/ or (3) reduced competition and predation in lower diversity rivers of the upland basins (Burns & Sidlauskas, 2019;López-Fernández et al., 2013). Furthermore, some of the large drainage basins adjacent to the modern Western Amazon, for example, the Upper Madeira and Paraguay basins, share many similarities in overall climate, water chemistry, and habitat physiognomy that may have facilitated geographic range expansions and population establishment .

| Species accumulation in lowlands
Initial diversification of curimatid lineages and phenotypes occurred in lowland Neotropical river basins during the Late Cretaceous, with diversification into clades with modern phenotypes and geographic ranges during the Paleogene and Neogene. This partitioning of the continent into core lowlands and peripheral uplands served to amplify the initial lineage diversity, generating the extraordinary levels observed in modern faunas. Low extinction rates have also been posited as a driver of high net diversification in NFF taxa (Lundberg et al., 1998), with most known paleofaunas having achieved modern phenotypes by the Late Miocene (Albert, Val, et al., 2018;Hoorn et al., 2010). However, extinction rate estimates are poorly constrained in the absence of a robust paleontological record (Beaulieu & O'Meara, 2015;Rabosky, 2010Rabosky, , 2016, and tropical rivers are poor environments for fossil formation and discovery (López-Fernández & Albert, 2011). Our hypothesis that Neotropical lowlands served as macroevolutionary museums is mainly supported by neontological evidence from the time-calibrated molecular phylogeny, showing accumulation of curimatid species in the lowlands (Figures 3-4).
Previous studies using macroevolutionary models have shown similar results. A study with the New World avifauna shows that tropical regions act as museums of diversity accumulating species diversity in low extinction rates (Gaston & Blackburn, 1996). In a molecular phylogenetic study of leaf beetles (Coleoptera: Chrysomelidae), McKenna and Farrell (2006) presented evidence that Neotropical forests are both museums and cradles of diversity showing distinct events of species diversification during the Cenozoic. The same pattern appears in other clades, as, for example, ants (Hymenoptera: Formicidae) that diversified during the Cretaceous were affected by the Cretaceous-Paleogene extinction, and experienced high diversification rates more during the Paleogene and Neogene (Moreau & Bell, 2013). More recently, a macroevolutionary study indicated bursts of speciation in three characiform families (Anostomidae, Serrasalmidae, and Characidae) around the Oligocene (~30 Ma), suggesting a cradle model of higher diversification rates driving the appearance of a remarkable species diversity in characoid fishes . The results found herein agree with those studies and emphasize the key role of the Amazon lowlands as a primary source of Neotropical biodiversity for many biotas (Antonelli, Ariza, et al., 2018;.

| Multiple dispersals to upland basins
Contrasting with lowlands, the cratonic uplands of the Guiana and Brazilian shields are better described as sinks of diversity, receiving multiple dispersal events of curimatid fishes during the Cenozoic, especially those of the most species-rich genera Cyphocharax and Steindachnerina, but also relevant dispersals of Curimata, Curimatella, and Pseudocurimata with closely related species in trans-Andean regions ( Figure 4). These multiple dispersals from lowland to upland basins occurred mostly during the Eocene and Oligocene, with a few events during the Neogene (Figures 3-4) associated with higher net diversification rates of widespread taxa ( Figure 5). These results reject previous hypotheses of origin and expansion of faunas from cratonic upland to lowland basins Eigenmann & Allen, 1942). Instead, our results support a pattern of biotic turnover potentialized by mega river captures in lowland basins that catalyzed subsequent events of allopatric speciation in large river basins draining the upland shields (Albert, Craig, et al., 2018;Lundberg et al., 1998).
The lowland origin in the POA followed by multiple independent colonizations of upland rivers is also supported by recent biogeographic reconstructions in other NFF clades. Long-whiskered catfishes (Pimelodidae), for example, emerged during the Late Cretaceous in the sub-Andean foreland basin, and had diversification directly affected by three major river captures between the Amazon and La Plata lowlands (Tagliacollo et al., 2015). Armored F I G U R E 5 GeoHiSSE range-dependent diversification plot of net diversification rates (M4, range-dependent process with hidden states). Colors of the branch outlines (blue to red) indicate the distribution of net diversification rates along the topology (<0. 001-0.26 (Fontenelle et al., 2021;Machado et al., 2018;Ramirez et al., 2020). In curimatids, lowland river captures between the Amazon/Paraguay and the Amazon/Tocantins (Figure 3, Figure   S1) facilitated the dispersal of taxa especially into the Brazilian Shield, although evidence for biotic or abiotic factors driving the upland colonization is not yet available, nor why those ancient lineages dispersed to uplands only in the Cenozoic.

| Comparisons with previous reconstructions and conclusions
Among several conclusions advanced by Vari (1988) bioregions. In addition, the results show, consistent with the vicariance model (Vari, 1988), that curimatids diversified in lowlands of the POA and Proto-Paraguay river basins, especially during the Early and Middle Miocene in the Pebas mega-wetland system of the modern Western Amazon (Albert, Val, et al., 2018;Hoorn et al., 2010). The role of lowland river captures between the POA and Proto-Paraguay leading to vicariance and geodispersal reinforces recent biogeographic reconstructions and interspecific relationships of NFF taxa from the Madeira and Paraguay basins Dorini et al., 2020;Hubert et al., 2007;Mateussi et al., 2017Mateussi et al., , 2020Melo et al., 2016;Ramirez et al., 2020;Tagliacollo et al., 2015), and the recent biogeographic reconstruction of more than 4760 Neotropical fishes despite small-scale differences . The discrepancy relative to this specific study is that curimatids are more informative on a larger biogeographic scale, likely due to the fewer narrowly distributed species compared with other species-rich groups such as characids (Characiformes), loricariids (Siluriformes), and rivulids (Cyprinodontiformes) .
Although the phylogenetic hypothesis of the present study is based on genetic information from six genes and about two-thirds of all known curimatid species, it is possible to test these results with studies using more complete taxon and gene sampling. For example, the inclusion of more lowland species of Curimatopsis might strengthen the signal of the lowland origin for the MRCA of all curimatids in the ancestral state reconstruction, the addition of Steindachnerina atratoensis will likely recognize the fifth event of dispersal from POA to the trans-Andean region, and the improvement of the phylogenetic resolution at the base of Cyphocharax sensu lato may refine the evolutionary and biogeographic history of that clade.
However, the addition of taxon/gene sampling is not expected to change the general results regarding lowland origins with multiple dispersals to uplands. Finally, our results demonstrate that molecular phylogenetic and parametric biogeographic studies can recognize similar biogeographical patterns when compared to more descriptive biogeographic studies (Vari, 1988). The coherence also reflects the predictive power of historical biogeography, which, regardless of the method, still exposes the intricate correlations among geomorphological processes and evolutionary relationships among Neotropical freshwater fishes.

ACK N OWLED G M ENTS
We thank the following people for thoughtful comments that con-

CO N FLI C T O F I NTE R E S T
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

DATA AVA I L A B I L I T Y S TAT E M E N T
Supporting Information is available online. Matrices, scripts, tables, trees, and other files are deposited in the Dryad Digital Repository: https://doi.org/10.5061/dryad.9p8cz 8wfw.