We used a global dataset of fish occurrence in discrete drainage basins to obtain species richness metrics (Leprieur et al., 2017; Tedesco et al., 2017a). These species richness values are point-estimates for fish “species inhabiting permanently or occasionally freshwater systems.” This dataset is subject to expected errors associated with biodiversity sampling and taxonomic knowledge. Given the size of the dataset, we do not expect these potential errors to bias the results.

Using shapefiles in Figure 2A we obtained basin-wide average topographic and climatic metrics (see Supplementary Dataset 1). For topography, we used the 90 m resolution, Shuttle Radar Topographic Mission (SRTM) digital elevation model downloaded from OpenTopography (Farr et al., 2007) to compute basin-wide average topographic metrics (see section “Basin-Wide Topographic Metrics”). For climatic metrics, we used the monthly mean, 30-year reanalysis dataset for the period of 1961–1990 with 0.5° spatial resolution (New et al., 2002) and the 29-year monthly mean air temperature data (UDel_AirT_Precip v4.01 product) for the period of 1981–2010 with 0.5° spatial resolution (Willmott and Matsuura, 2001) provided by the Physical Sciences Laboratory (National Oceanic and Atmospheric Administration, NOAA).

Given that tectonic activity impacts biodiversity, we created two data partitions such as tectonically active or tectonically stable, which we identified based on seismic activity. Seismic activity was identified through Peak Ground Acceleration (PGA) data obtained from the Global Earthquake Model (GEM) (Pagani et al., 2018). GEM uses several methods to obtain PGA, including compilation of hazard maps and reduction of seismic data. In this dataset, PGA is the 10% probability of exceeding the 50-year reference shear wave velocities (see Pagani et al., 2018). We use a threshold of 3.2 to distinguish between tectonically active vs. stable regions based on the observation that the Amazon basin contains an average PGA of 3.2 (no-data grid cells are not averaged) and is a tectonically stable region.

Tedesco et al. (2017a) provides a list of fish species occurring in freshwater for 3,119 basins globally out of which we were able to acquire topographic data for 3,038. For every basin in the dataset, we calculated species richness (SR) as the total number of valid fish species (Figures 2C,D). Species density (SD) was then calculated via a regression of SR on drainage area. The Species-Area Relationship (SAR) assumes the form: SR = Area^z (MacArthur and Wilson, 1966) and, therefore, SD was calculated as SR/Area^z, where the exponent z is obtained from the power-function regression of SR on Area (Rosenzweig, 2004; Albert et al., 2011). We then created data partitions based on seismic activity (active and stable) and latitude (tropics and extra-tropics) to assess differences in SAR based on geologic and climatic settings. Lastly, we focused on tropical regions as basins between 23.5 degrees latitude north and south, where most fish species live. We examined the effects of topography by separating basins into uplands and lowlands defined as median elevation above and below 500 m above sea level, respectively. We then assessed the importance of climatic and topographic predictor variables such as precipitation and relief, respectively, in individually predicting SD in the tropics. We create four sub-groups: (1) tectonically active highlands; (2) tectonically stable highlands; (3) tectonically active lowlands; (4) tectonically stable lowlands.

We used TopoToolbox to extract topographic metrics for each drainage basin (Schwanghart and Scherler, 2014). We obtained the average values for elevation above sea level as well as topographic relief. The latter was computed as the absolute range of topography over a 2,500 m moving window at every cell within a given basin.

Relief is a topographic metric that describes the local amplitude of topography. In tectonically active regions, relief scales directly with uplift rates and is often used as a proxy to identify relative differences in tectonic uplift (i.e., the rate of advective motion of rock) (e.g., Montgomery and Brandon, 2002; Kirby and Whipple, 2012). Relief also scales with a rock’s resistance to erosion; harder (or softer) rocks promote steeper (or gentler) rivers (e.g., Hack, 1973; Duvall et al., 2004; Gallen, 2018). Thus, actively uplifting/eroding regions with highly variable lithology promote complex transient evolution and topography (e.g., Forte et al., 2016). Importantly, landscapes that are actively changing due to some past perturbation (i.e., base-level fall, river capture, or tectonic uplift) will have a positive correlation between relief, river steepness, and erosion rates within a given basin and, therefore, the amount of sediments actively fluxed through rivers (Portenga and Bierman, 2011; Kirby and Whipple, 2012; Gallen et al., 2013). Based on a previous modeling study, the ratio of the magnitude of the perturbation to the initial landscape relief dictates the degree of drainage reorganization, which in turn affects the placement or removal of dispersal barriers for aquatic organisms (Lyons et al., 2020; Stokes and Perron, 2020). Thus, average basin relief is a good metric for overall topographic steepness of a river basin and for linking topographic responses to autogenic and exogenic forcings. Also, it has been directly assessed with riverine species evolution in modeling studies (e.g., Lyons et al., 2020). Nonetheless, we emphasize relief is the time-integrated outcome of climate, tectonics, and surface and groundwater processes, and not strictly the outcome of landscape transience.

The SD values of freshwater fishes among river basins of the world vary systematically by seismic setting and latitude. We found a significant relationship (R² = 0.79; p < 0.01) between species richness (SR) and species density (SD) among all river basins worldwide. We found a relatively weak although significant relationship (R² = 0.29; p < 0.01) between SD and area among the largest basins (> 100,000 km²; n = 133), but not for the medium (> 10,000–100,000 km²; n = 504) or smaller basins (< 10,000 km², n = 2,482).

For both tropical and extra-tropical regions we observe that tectonically stable settings have a higher baseline species density value (i.e., y-intercepts of the SAR) but a lower slope (Figures 3, 4). Moreover, the set of tropical basins exhibits a higher slope value than does the set of extra-tropical basins. Estimated this way, there is no significant relationship between fish species density and drainage area for all of the world’s freshwater basins, although there is a modest relationship between these variables among the 133 basins larger than 10,000 km² (R² = 0.12; p < 0.01).

The most species-dense drainage basins are substantial outliers in each of the geographic sets of basins by tectonic activity and latitude (Figures 3, 4). The statistical distributions of species density within each set reinforce the differences between tectonically stable and active regions (Figure 5). In all cases, the distributions are approximately log-normal, with the highest species densities with heavier tails in the tectonically stable group (Figures 5A,B). The average species density is higher in the tectonically stable regions compared to active areas irrespective of latitude (Figures 5C,D).

Elevation (i.e., uplands and lowlands) is a well-known parameter controlling the distribution of aquatic species (Albert et al., 2018a). As expected, lowlands contain a 10-fold higher baseline species density (i.e., intercept of the SAR) but a lower slope. Basins affected by Neogene mega river captures (i.e., Amazon, Congo, Mekong) are exceptionally diverse and outliers in the SAR (Figures 2B, 6).

In this study, we focus on river basins situated in the tropics from where a majority (66%) of freshwater fish species are known. Based on the observation that uplands and lowlands have differing relationships with topographic metrics depending on tectonic activity, we grouped basins based on drainage area, mean elevation, and tectonic activity and assessed four groups: tectonically active and stable lowlands and uplands.

Absent rainfall rates, we find no systematic relationships between species density and temperature, elevation, and relief when considering all basin sizes in each data subcategory. However, we find statistically significant correlations with these metrics for basins greater than 10,000 km² which is an approximate threshold at which landscape evolution processes might impact biodiversity (Albert et al., 2018a,2021). Given this observation, we describe the following results for basins larger than this drainage area threshold.

The results of regressions of freshwater fish SD values against two climatic (i.e., precipitation and temperature) metrics among river basins of the world are consistent with general macroecological expectations of highest diversity in wetter, warmer regions (Worm and Tittensor, 2018). Among freshwater basins globally, a disproportionate number of fish species inhabit the lowlands of large tropical river basins in South America, Africa, and Southeast Asia, under predictable ecological conditions of large geographic area, warm, humid tropical climate, relatively flat topography (i.e., low average relief), and high habitat volume (i.e., high precipitation and run-off) (Lundberg et al., 2000; Oberdorff et al., 2011; McGarvey and Terra, 2016; Leprieur et al., 2017; Antonelli et al., 2018b). The six river basins with greatest fish SR values (i.e., Amazon, Congo, La Plata, Orinoco, Mekong, and Tocantins basins) have a combined total of more than 9,200 species in an area of 13.7 million km², thus accounting for 51% of all freshwater fishes globally, in an area of just 14% of the 98.7 million km² of all river basins on Earth combined.

Among both tropical and extra-tropical basins, the species-area exponent is higher for tectonically (seismically) active than stable regions (Figures 3, 4). This result suggests higher rates of dispersal among basins of low-relief tectonically passive margins, as indicated by previous studies showing higher SR in tectonically stable settings (Badgley et al., 2017; Griffiths, 2018). The higher correlation between area and SD for stable basins in tropics than extra-tropics indicates more heterogeneous geographical conditions and geological history of the many (n = 1,478 or 47% of the total) river basins in extra-tropical and tectonically stable regions (Figure 4).

The observed negative relationship between SD and average basin elevation among tectonically active uplands basins is consistent with previous observations of species occurrence along elevation gradients in the Amazon (Lujan et al., 2013), and across the South and North American continents (Smith et al., 2010Griffiths, 2018). Relationships of SD with air temperature and precipitation in these regions (Figures 7A,B) are consistent with contributing roles of both contemporary ecological and historical (time-integrated) macroevolutionary effects driving down aquatic species diversity at higher elevations (e.g., Lujan et al., 2013; Hazzi et al., 2018). Similarly, the history of topographic growth, as opposed to present-day topographic or climatic conditions, is thought to be more important to the evolution and enrichment of many terrestrial (Castroviejo-Fisher et al., 2014; Antonelli et al., 2018a; Azevedo et al., 2020; Réjaud et al., 2020; Igea and Tanentzap, 2021) and aquatic (Smith et al., 2010; Badgley et al., 2017) vertebrate faunas.

Macroevolutionary theory predicts that regional SR values arise from interactions among three fundamental parameters: rate of net diversification (δ), time of net diversification (t), and regional carrying capacity (S_max) (Rabosky, 2010). The per-species net diversification rate (δ) is: δ = λ + d—μ, where per-species speciation (λ) and dispersal (d) rates add species, and per-species extinction rate (μ) removes species from a region. When diversity is unbounded (i.e., δ is independent of S_max), δ and t provide limits to diversity; i.e., δ = δ0, where δ0 is the intrinsic diversification rate. Under these non-equilibrium conditions, the number of species at time t is: SRt = e^δ0⋅t (Cornell, 2013). Such a model applies when speciation and dispersal rates are low relative to total available niche space and/or geological age of a region, or when μ > λ + d for sufficient time that a diversity limit (i.e., SR_max) is not approached. When diversity is bounded, SRt depends on SR_max and t; i.e., δ = δ0—aSR, where a is the strength of diversity-dependent feedback on δ, the carrying capacity (S_max = δ0/a) depends on the time-integrated δ or δt = ∫ [λ(t) + d(t) – μ(t)] dt, and expected equilibrium species richness is: SRt = e^{δ t} (Rabosky, 2013; Cornell and Harrison, 2014; Harmon and Harrison, 2015; Rabosky and Hurlbert, 2015).

Under this theory, the smaller difference in correlation values between tectonically active basins (5%) among tropical vs. extra-tropical (9%) basins is unexpected (Stanley, 2014; Albert et al., 2017), because there are many more tectonically stable basins (Figure 4). This result is partly due to the much wider range of SD values among large basins on stable terrains, in particular from the many large basins from cold boreal regions (e.g., northern Canada and Russia) and arid tropical regions (e.g., northern Africa, central Australia) with low SD values (Figure 1 lower panels). There are fewer counterparts of these large low SD basins on active terrains.

Both macroecological and macroevolutionary models predict higher SD values in the dendritic geometry of river drainage networks (Fagan, 2002; Thomaz et al., 2016). River networks are a more spatially fragmented substrate than an equivalent Euclidean landscape with the same surface area (Rodriguez-Iturbe and Rinaldo, 2001; Dias et al., 2013). The hierarchical-branching of drainage networks is a more effective geometry for breaking up a geographically widespread species into daughter species (Wiens, 2002; Muneepeerakul et al., 2007; Tonkin et al., 2018). However, even more effective than a 3D dendritic surface for fragmenting and merging populations is a 4D dendritic surface changing in time; i.e., river capture (Albert and Crampton, 2010; Albert et al., 2017, 2018a).

The results of this study are consistent with several predictions particular of the RCH (Figure 11; Albert et al., 2018a). By merging geographic areas (geodispersal), river capture facilitates organismal dispersal and gene flow, and therefore acts to slow rates of speciation and extinction, i.e., lower species turnover (Albert and Crampton, 2010). However, by subdividing areas (vicariance) river capture also acts to increase rates of speciation and extinction. The results of this study suggest that genial ecological conditions are necessary, but not fully sufficient, to explain the basins with the highest diversity, consistent with the prediction of the RCH that dispersal across the watershed margins of adjacent lowland basins increases basin-wide SR and SD values. These increases occur both by adding individuals of different species, and adding individuals of existing species, thereby lowering the within-basin extinction rate (i.e., rescue effect; Brown and Kodric-Brown, 1977).

Dispersal of freshwater fishes among adjacent river basins may occur by multiple abiotic and biotic mechanisms (e.g., Tagliacollo et al., 2017). In some cases, the effects of geodispersal by river capture and biotic dispersal by organismal movements can be difficult to separate, for examples in seasonally flooded wetlands that straddle low-elevation drainage divides; e.g., Rupununi and Izozog swamps at the margins of the Amazon and adjacent basins which are sites of longer-term river captures and seasonal dispersal (Albert et al., 2011). In other cases, riverine corridors serve as ecological filters in which organismal trait values (adult body size, habitat utilization, tropic specializations) influence the species richness and composition of the biotic interchanges due to river capture; e.g., Casiquiare River.

The Mega Capture Hypothesis (MCH) predicts that large river captures (> 10,000 km²) leave a disproportionately enhanced signature on the accumulation of basin-wide SR values in riverine and riparian taxa (Tagliacollo et al., 2015; Albert et al., 2018a,2021). The basins with highest fish SD values (i.e., Amazon, Orinoco, La Plata, Congo, Mekong) have all been exposed to the effects of mega river captures within the past 20 Ma, which merged portions of the whole biotas of riverine and riparian taxa among adjacent basins (e.g., Bragança and Costa, 2019; Musher et al., 2019; Van Steenberge et al., 2020; Albert et al., 2021; Chen et al., 2021; Sun et al., 2021; van der Merwe et al., 2021).

Although SR and SD are significantly correlated among basins worldwide (R² = 0.8, n = 3,038, p < 0.01), basins with highest SD values drain primarily tectonically stable regions and these basins were assembled during the Neogene (c. 23–2.6 Ma) and Quaternary (2.6–0 Ma) through the action of mega river capture events (Albert et al., 2018a). Among the largest basins worldwide (i.e., those > 100,000 km²), basins with the top 10 SD values are (in descending order): Amazon, Orinoco, Chao Phraya (Thailand), Mekong, Essequibo, Paraná-Paraguay, Congo, Tocantins, Uruguay, and Zhujiang (Pearl) basins. Under the RCH, the unexpectedly high SD values of these river basins, as assessed by their positive deviations from the regression in Figures 3–6, arose from the merging of multiple older and smaller basins through mega river capture events (Albert et al., 2018b,2021; Sun et al., 2021). All these basins have undergone substantial changes to their watershed margins over the last 20 million years, with significant portions, sometimes>50% (e.g., Hoorn et al., 2010) of their modern areas having been affected by river capture during this time (e.g., Clark et al., 2004; Goudie, 2005; Hoorn et al., 2010).

The Intermediate Capture Rate Hypothesis (ICH) predicts highest SD values on landscape with “just right” rates of river captures through time and at appropriate spatial scales (Albert et al., 2018a). The results of SD regressions against two landscape metrics (i.e., elevation and relief) are consistent with several predictions of the ICH. Maximum SD values are obtained in tectonically stable lowland basins (Figure 10), where river capture dynamics are expected to drive an excess of speciation and dispersal events as compared with extinction events per unit time, and therefore a net accumulation of SR through time. By contrast, SD values are lower, and not correlated with relief, among seismically active lowland basins (Figure 9) and stable upland basins (Figure 8). Finally, SD values are lowest, and are negatively correlated with relief, among the set of active uplands basins (Figure 7).

Basins with highest fish SD values are located within the tropics and on stable terrains, although not all basins in the tropics or on stable terrains have high fish species density (Figure 2). This is similar to the results of Albert et al. (2018a) who show differing river capture rates on stable vs. active tectonic platforms. Rates of river capture are poorly constrained on most landscapes worldwide, but preliminary evidence from South America indicates they may be several orders of magnitude faster on alluvial lowland sedimentary basins of continental interior than on upland cratonic regions (Ruokolainen et al., 2019; Goldberg et al., 2021).

Landscape evolution modeling results arrive at similar conclusions and suggest that river captures are likely more frequent and larger in low-relief landscapes (Lyons et al., 2020), which exhibit highest SD values worldwide (Figure 2B). Conversely, high relief in tectonically active landscapes acts to fragment the species range and increase extinction rates at smaller spatial scales as compared to low-relief landscapes (e.g., Smith et al., 2002; Albert et al., 2006, 2018a; Borregaard et al., 2012; Griffiths, 2018). Binned by relief, our dataset supports these directional relationships and reveals that SD values peak at an intermediate relief value when comparing tectonically stable lowlands with tectonically active uplands (Figure 12). The probability of river captures increases where neighboring basins erode laterally at different rates, which especially true near topographic escarpments with uniform rock types (Salgado et al., 2014; Willett et al., 2018; Calegari et al., 2021; Wang and Willett, 2021). The across-divide differences in mean relief, elevation, slope, and other topographic metrics that measure steepness, dictate the direction and rate at which divides migrate (Whipple et al., 2017). Importantly, in lower relief settings, landscape perturbations by long-wavelength (100 s of km) and low-amplitude (<1 km) uplift (i.e., dynamic topography; e.g., Bicudo et al., 2019), local uplift (i.e., faults), and base-level fall are likely to reach or surpass the observed ranges of relief, which more easily prompts drainage reorganization (Lyons et al., 2020). To the extent that the value of topographic relief affects the frequency of river captures across these landscape settings (e.g., Lyons et al., 2020), our findings are consistent with the expected effects of barrier displacement caused by landscape transience and support the ICH (Albert et al., 2017).

According to the macroevolutionary model outlined above, we may expect to see SR_max values at intermediate rates of river capture under non-equilibrium conditions, when basin-wide species richness (SR) values are growing because the rate of speciation (λ) exceeds extinction (μ). However, under more equilibrium conditions when the rates of speciation (λ) and extinction (μ) are similar, SR may be more strongly controlled by the regional carrying capacity (SR_max; Rosenzweig, 2004; Albert et al., 2017). Results of this study indicate that rates of river capture on stable lowland continental platforms are sufficiently slow enough to allow speciation to occur among isolated portions of river networks through time (t sensu Figure 1), while also being sufficiently fast enough to allow dispersal to populate adjacent basins. This combination of parameter values allows SR values to increase in lowland sedimentary basins which experience higher rates of river capture and reduce the extinction risk of already resident species (Fagan, 2002; Tedesco et al., 2012; Douglas et al., 2013).

By contrast, slower river capture rates on tectonically stable and more erosion resistant upland regions (e.g., continental cratons) are expected to inhibit dispersal among adjacent basins, and thereby lower the rate of increase of diversity through time (lower δ sensu Figure 1). By the same logic, faster rates of river capture within alluvial lowland sedimentary basins of continental interior or coastal plains retard fish diversification, because populations do not become isolated for long enough to allow genetic divergence. In positing that lineage diversification dynamics arises from the multiple effects of dispersal and gene flow on speciation and extinction, the ICH resembles the Shifting Balance Theory (Wright, 1982) and Effect Hypothesis (Vrba, 1983), with the notable differences that, under the ICH, speciation may occur due to genetic drift alone rather than requiring natural selection, and species may therefore not necessarily be adapted to different ecological niches (Harvey et al., 2019). The hypothesis that highest SR values are obtained at intermediate capture rates reflects a larger perspective that all possible evolutionary drivers impose trade-offs on organismal diversification, achieving maximal effectiveness over a limited domain of parameter values; e.g., the intermediate disturbance and productivity hypotheses (Huston, 1994; Fraser et al., 2015).

The log-normal SD frequency distributions observed in this study often characterize biodiversity profiles and other biological systems that grow over time from the multiplicative interactions of many independent random variables (Crow and Shimizu, 1987; Rozenfeld et al., 2008; Magurran, 2013). As numerous studies suggest that landscape transience, river captures, and escarpment migration are common characteristics of intracontinental lowland regions (e.g., Harbor et al., 2005; Gallen et al., 2013; Val et al., 2014; Beeson et al., 2017; Gallen, 2018; Willett et al., 2018; Wang and Willett, 2021), landscape evolution processes might be a common underlying mechanism of diversification in continental regions.

This study examines relationships of fish species richness with possible drivers among basins assigned to broad latitudinal categories (i.e., tropical and extratropical), but does not examine possible effects of latitude on habitat heterogeneity within river basins. Such an analysis would be complicated by many additional factors, with possible expectations for greater habitat heterogeneity in tropical than extra-tropical basins, for basins with N-S than W-E main-stem axis orientations, for upland than lowland basins, and for stable than active terrains. Many other potentially important factors could also be examined, including especially distance from continental geographical centroid or center of connectivity (Smith et al., 2010), and mean or maximum phylogenetic clade age (Miller and Román-Palacios, 2021).

This study also uses topographic metrics such as relief as proxy for landscape transience. Spatial variability in relief is not a unique outcome of exclusively landscape transience. The erosive susceptibility inherent to lithologic types, for example, also influences relief. The primary control of relief and its strong correlation with erosion has been recognized since the early days of geomorphology (Gilbert, 1877), and this correlation continues to be identified using state of the science techniques to measure erosion rates (e.g., von Blanckenburg, 2005). The cross-divide difference in erosion rates can be especially indicative of transient river network reorganization (e.g., Willett et al., 2014; Whipple et al., 2017). Future studies can incorporate cross-divide relief and erosion differences with computational tools such as Forte and Whipple (2018) along with our approach to further investigate links among river captures and species richness.