Threats of global warming to the world’s freshwater fishes

1Freshwater habitats are disproportionally biodiverse. While they cover only 0.8% of the Earth’s surface, they host ~15,000 fish species, corresponding to approximately half of the global known fish diversity,. Freshwater habitats are also disproportionally threatened by human activities and environmental change, which have resulted in substantial declines in freshwater biodiversity over the past decades,. Amid human pressures on freshwater ecosystems (including water abstraction, diversion, damming, and pollution), anthropogenic climate change is expected to become increasingly important in the future,. Rising air temperatures and changing precipitation patterns modify water temperature and flow regimes worldwide, thus affecting two key habitat factors for freshwater species. Being ectotherms, fish are directly influenced by water temperature, while the hydrologic regime determines the structure and dynamics of the freshwater habitat,. In addition, the insular nature of many freshwater habitats may hamper compensatory movements to cooler locations, especially for fully aquatic organisms like fish. Recent continental and global studies have underscored the high vulnerability of freshwater fish species to climate change–. Yet, potential impacts of climate change on freshwater fishes have not yet been comprehensively assessed, in sharp contrast with the many studies assessing potential climate change impacts on species in terrestrial systems–.

2Here, we assess future climate threats to 11,425 riverine fish species by quantifying their exposure to flow and water temperature extremes under different global warming scenarios. We focus on extremes rather than hydrothermal niche characteristics in general, because extremes are more decisive for local extinctions and potential geographic range contractions,. Following the latest IPCC report, we include scenarios that limit global mean temperature increases to 1.5 and 2.0 °C. For comparison purposes, we include two additional scenarios: a “current pledges” scenario set at 3.2 °C warming and a “no-policy” scenario (no mitigation) set at 4.5 °C warming (all temperatures relative to pre-industrial),. The 3.2 °C warming scenario represents the maximum warming predicted to occur by the end of the century (with 66% probability) if all current greenhouse gas emissions reductions targets (unconditional Nationally Determined Contributions) for 2030 are met and no further cuts are performed. We calculate the present and future weekly flow and water temperature values corresponding with each warming level at a spatial resolution of 5 arcminutes (~10 km) using a global hydrological model coupled to a dynamic water temperature model,. We force the hydrological model with meteorological input from five Global Climate Models (GCMs) combined with four Representative Concentration Pathway (RCP) representing future greenhouse gas emissions. To assess the threat imposed by future climate extremes, we first retrieve flow and water temperature extremes corresponding with the current climate. Across the geographic range of each species, we quantify the maximum and minimum weekly water temperature and flow, as well as the number of zero flow weeks (see “Methods”). We then define the magnitude of threat for each fish species in a given global warming scenario as the proportion of the geographic range of the species where projected extremes in water flow and temperature will exceed those defined based on the current climate within its ranges. We do this for two dispersal assumptions: ‘no dispersal’ assuming that each fish species is restricted to its current geographical range, and “maximal dispersal” assuming that each fish species can move beyond its current range within a surrounding region delineated by the intersection of watersheds (hard physical boundaries) and freshwater ecoregions (i.e., regions with similar evolutionary history and other ecological factors relevant to freshwater fishes). Finally, we use phylogenetic regression relating the species’ threat levels (i.e., the proportion of the range exposed to future climate extremes beyond current levels within the range) for each warming level and dispersal assumption to a suite of relevant species characteristics.

3We find clear differences in the magnitude of threat between the different warming scenarios. In a 3.2 °C warmer world, 36% of the species have over half of their present-day geographic range exposed to climatic extremes beyond current levels (no dispersal assumption). This number reduces to 9% of the species in a 2 °C warmer world and to 4% of the species if warming is limited to 1.5 °C. We conclude that for protecting freshwater biodiversity, commitments to limit global warming need to be strengthened.

4The scenario without climate-change mitigation policy (+4.5 °C) and without dispersal resulted in at least half of the geographic range threatened by projected climate extremes for 63% (±7%) of the freshwater fish species. Assuming maximal dispersal for the same warming level, the proportion of species with over half of their geographic range threatened decreased to 24% (±13%). The values in brackets represent the standard deviation of the GCM–RCP combinations ensemble for that warming level and dispersal assumption (Supplementary Fig. 1). The proportion of species with more than half of their range threatened was projected to decrease to 8–36% (±3–11%), 1–9% (±1–4%), and 1–4% (±0–2%) for warming levels of 3.2 °C, 2 °C and 1.5 °C, respectively, with the larger values for the no dispersal assumption (Fig. 1).

5We found hotspots of future climate threat in tropical, sub-arid and Mediterranean regions (Fig. 2). At low warming levels, hotspots are restricted to small areas within tropical South America, North-East Mexico, southern US, southern Europe, Southern Sahara, central Africa (large lakes), Middle-East, India–Pakistan, South-East Asia, and western Australia. At higher warming levels, hotspots are considerably larger, particularly in South America, southern Europe, India–Pakistan, and Australia. At higher latitudes, threats become prominent only at higher warming levels (3.2, 4.5 °C). Overall, threats are largest in tropical watersheds such as the Amazon, Parana, Tocantis, Niger, Senegal, Zambezi, and Chao-Phraya (Fig. 3; see Supplementary Fig. 4 for a more exhaustive overview). Watersheds in non-tropical areas characterized by relatively high threat levels are the Don and the Danube in Europe, and several watersheds in Australia (Fig. 3). Under the maximal dispersal assumption, locations of threatened areas are similar to those under the no dispersal assumption but with lower threat levels than in the no dispersal assumption (Figs. 2 and 3).

6Our findings indicate that freshwater fish species are primarily threatened by climate change-induced increases in maximum water temperature, whereas amplified extreme flow conditions are considerably less important (Fig. 4; Supplementary Fig. 2 and Supplementary Tables 3, 4). Projected reductions in minimum water temperature pose virtually no threat (Supplementary Fig. 3 and Supplementary Tables 3, 4). Therefore, the spatial patterns of future climate threat mostly resemble the patterns of threats due to increased maximum water temperature (Fig. 4). Areas affected by changes in low flow are mainly observed in upstream reaches in South America, parts of the central US, around the Mediterranean Sea and Middle East (Fig. 4 and Supplementary Fig. 3). This reflects that climate change is projected to result in more severe low flow conditions mainly in drought-prone regions, while water temperature rises almost everywhere (Supplementary Fig. 5-V). In addition, changes in low flow conditions might be more relevant for smaller upstream streams not captured within the ~10 km grid-cell resolution of the global hydrological model employed in this study. In contrast, areas affected by changes in high flow are confined to a few downstream segments of the main stems of large rivers (Supplementary Fig. 3 and raster layers provided as Supplementary Data 2 and 3). Our results further show only limited overlap of threats imposed by amplified flow and water temperature extremes, reflecting the dissimilar spatial distribution of both threats (Fig. 4 and Supplementary Tables 3, 4).

7According to our phylogenetic regression models (n = 9,779 species), the magnitude of threat imposed by future climate extremes is mainly related to species’ habitat type and current geographic range size, followed by IUCN threat status and body length (Fig. 5). Threats are much lower for species that live across the freshwater and marine realms (note that our projections concern the freshwater environment only, thus ignoring potential climate threats within the ocean; Supplementary Table 5). In line with this, relatively small threat levels are found for orders mostly comprising diadromous species, such as Mugiliformes (mullets), Osmeriformes (smelts), Syngnathiformes (e.g., pipefish), Tetraodontiformes (e.g., pufferfish) and Pleuronectiformes (flatfish) (Supplementary Fig. 6). Further, species with a smaller geographic range and body size are more likely to be threatened by climate change (negative regression coefficients; Supplementary Table 5). We also noticed lower threat levels for species currently belonging to a low IUCN threat category (e.g., “near threatened” or “least concern”; Supplementary Table 5). We found similarly low threat levels for species that are “data deficient” within the IUCN Red List (43% of the 9,779 species analyzed). Future climate threats were only weakly related to climate zone, commercial importance category and trophic group (Fig. 5). The results of the traits analysis were largely consistent between the two dispersal scenarios (Fig. 5). The results were less consistent across the warming levels, whereby the importance of geographic range size dropped considerably at higher warming levels under the no dispersal assumption, while habitat type became more important (Fig. 5).

8This study represents the first comprehensive assessment of the threat of potential future climate extremes to freshwater fish species, covering both flow and water temperature, the entire globe and about 90% of the known freshwater fish species. We found that in a “current pledges” scenario (3.2 °C warming), over one third of the freshwater fish species is projected to have more than half of their geographic range threatened by future climate extremes beyond current levels. Threats are considerably reduced under the assumption of maximal dispersal, yet this might represent an overly optimistic estimate given current and future barriers to dispersal. Our results suggest that increases in maximum water temperature constitute a larger threat to freshwater fishes than changes in minimum water temperature or high and low flow conditions. This is because water temperatures vary less within species ranges and are projected to rise almost everywhere, while flow conditions are more spatially variable hence projected future flow is less likely to exceed present-day extremes within the species’ ranges. In line with previous studies, we found that climate change will result in reduced flows mainly in drought-prone regions,. In addition, depletion of low flows might be most important at low stream orders, which are not well captured by the 5 arcminutes resolution of the hydrological model PCR-GLOBWB employed in this study. While global estimates of hydrological variables are available at higher spatial resolutions, 5 arcminutes is the highest resolution currently achievable for future projections of both flow and water temperature. Hence, the spatial resolution of our analysis might result in an underestimation of the impacts of climate change on species living in smaller upstream reaches. We further note that we did not explicitly consider changes in flow or water temperature seasonality, which might disfavor species whose life histories are adapted to specific flow or temperature regimes (e.g., specific seasonal flow regimes),. However, flow and water temperature seasonality were clearly correlated to the minimum flow and water temperature extremes within the species’ ranges (Pearson’s r of 0.6–0.9), indicating that the effects of changes in the seasonality of flow or temperature were at least partially covered by our predictions.

9Exposure to climate extremes beyond present-day values does not necessarily imply local extinction. If species’ current distributions are confined by factors other than flow or water temperature (e.g., biogeographic dispersal barriers or anthropogenic pressures), species might be able to withstand larger temperature and flow extremes than inferred based on their current geographic range,. The same holds if species are able to adapt to new water temperature and flow conditions or if fishes have the possibility to hide from extremes in micro-climatic refugia, for example due to water stratification or small-scale thermal heterogeneity, which are not included in our hydrological model. On the other hand, species’ range maps are relatively coarse representations of species occurrence, hence some species might be more affected than indicated by our results (i.e., if present-day climate extremes within their geographic range already preclude local occurrence). Indeed, a tentative comparison of our species-specific maximum weekly water temperature limits with critical thermal maxima reported from laboratory tests suggest both under- and overestimations by our geographic range-based thermal limits, while showing an overall reasonable agreement (mean percent difference = 9%; Pearson’s r = 0.62; Supplementary Fig. 7). Further work is required to better understand deviations between our empirical thresholds and the lab-based maxima, which may stem not only from uncertainties in our modeling approach (e.g., in the range maps or the water temperature model) but also from heterogeneity in experimental conditions. Additionally, we recognize that a given increase in global mean temperature may lead to locally different exposure levels depending on the GCM–RCP combination, as each is characterized by specific distributions of changes in water temperature and flow. We notice a greater variability across GCMs than RCPs when looking at species-specific proportions of range affected (Supplementary Fig. 1), similar to previous findings for hydro-climatic variables. However, variability across the GCMs did not affect the species-specific thresholds, which were consistent across the models (Supplementary Fig. 8).

10Although our future climate threat assessment is associated with uncertainties, our comparative analysis across the different warming levels and targets clearly showed a sharp increase in potential impacts with increasing global warming. Species already listed as “endangered” or “critically endangered” on the IUCN Red List of threatened species might be particularly affected by future warming, as these species were characterized by the highest future climate threat levels. Our findings also show that threats imposed by amplified climate extremes are expected to be particularly high in tropical watersheds, in accordance with previous studies suggesting large climate-change induced freshwater habitat degradation in the tropics–. Tropical species are indeed expected to be highly affected by climate change, and our results confirm this for freshwater fishes. Many tropical watersheds host low-income food-deficit countries where local communities are highly dependent on fishery as a primary food source. Indeed, up to 50% of household incomes in countries along the Mekong, Zambezi, and Brazilian Amazon depend on fishing. Hence, increased exposure of freshwater fish species to climate extremes, potentially resulting in local extinctions, may have important socio-economic repercussions in these regions. Our findings indicate that limiting global warming to 2 °C will reduce the proportion of freshwater fish species with more than half of their range threatened by 74–81% (the range refers to the two dispersal scenarios) compared to current pledges of governments (3.2 °C). Restricting the global mean temperature rise to 1.5 °C will lower this proportion by an additional 11–14% (or 53-58% compared to 2 °C). While we acknowledge that the ecological realism of our model projections can be improved, these first comparative estimates highlight the need to intensify (inter)national commitments to limit global warming if potentially severe disruption of freshwater biodiversity is to be prevented.

11We compiled species’ geographic ranges from a combination of datasets. We employed the IUCN Red List of Threatened Species database, which provides geographic range polygons for 8,564 freshwater fish species (~56% of freshwater fish species), compiled from literature and expert knowledge. We complemented these ranges with data from Barbarossa et al, who compiled geographic range polygons for 6,213 freshwater fish species not yet represented in the IUCN dataset, and the Amazonfish dataset, which provides range maps for 2,406 species occurring in the Amazon basin. We harmonized the species names based on Fishbase (www.fishbase.org) and merged the ranges (i.e., union of polygons) from the different datasets to obtain one geographic range per species. We then resampled the range polygons of each species to the 5 arcminutes (~10 km) hydrography of the global hydrological model (see below), with a given species marked as occurring in a cell if ≥ 50% of the cell area overlapped with the species’ polygon. In total, we obtained geographic ranges for 12,934 freshwater fish species, covering ~90% of the known freshwater fish species. We excluded 1,160 exclusively lentic species because our hydrological model is less adequate for lakes than for rivers, i.e., it does not account for water temperature stratification (see section “Phylogenetic regression on species traits” for an explanation of how habitat information was extracted). Out of the 11,774 (partially or entirely) lotic fish species, we excluded 349 species (~3%) because their occurrence range was smaller than ~1,000 km2 (i.e., ten grid cells), which we considered too small relative to the spatial resolution of the hydrological model (see below). Hence, the analysis was based on 11,425 species in total (Supplementary Figs. 9, 10; a raster layer providing the number of species at each five arcminutes grid cell is available as Supplementary Data 6).

12We employed the Global Hydrological Model (GHM) PCR-GLOBWB with a full dynamical two-way coupling to the Dynamical Water temperate model (DynWAT) at 5 arcminutes spatial resolution (~10 km at the Equator), to retrieve weekly streamflow and water temperature worldwide,. PCR-GLOBWB simulates the vertical water balance between two soil layers and a groundwater layer, with up to four land cover types considered per grid cell. Surface runoff, interflow, and groundwater discharge are routed along the river network using the kinematic wave approximation of the Saint–Venant Equations and includes floodplain inundation. Apart from the larger lakes, PCR-GLOBWB includes over 6,000 man-made reservoirs as well as the effects of water use for irrigation, livestock, domestic, and industrial sectors. PCR-GLOBWB computes river discharge, river and lake water levels, surface water levels and runoff fluxes (surface runoff, interflow and groundwater discharge). These fluxes are dynamically coupled to DynWAT along with the meteorological forcing, such as air temperature and radiation from the GCMs to compute water temperature. DynWAT thus includes temperature advection, radiation and sensible heating but also ice formation and breakup, thermal mixing and stratification in large water bodies, effects of water abstraction and reservoir operations. We selected this model combination because it allows a full representation of the hydrological cycle (considering also anthropogenic stressors, e.g., water use), it fully integrates water temperature and calculates the hydrological variables on a high-resolution hydrography. The choice of one hydrological model over an ensemble was motivated by the fact that very few GHMs or Land Surface Models calculate water temperature at the spatial resolution desired for this study,. The PCR-GLOBWB model setup was similar to Wanders et al., with the exception that flow and water temperature were aggregated at the weekly scale to capture the fish species’ tolerance levels to extreme events.

13To assess climate change threats to freshwater fishes, we focused on climate extremes rather than hydrothermal niche characteristics in general, because extremes are more decisive for local extinctions and potential geographic range contractions,. We quantified climate extremes using long-term average maximum and minimum water temperature (Tmax, Tmin), maximum and minimum flow (Qmax, Qmin), and the number of zero flow weeks (Qzf), based on the weekly hydrograph and thermograph of the hydrological model. Water temperature is considered the most important physiological threshold for fish species, as mortality of ectothermic species occurs above and below lethal thresholds,. Decreases in minimum flow directly affect riffle-pool systems and connectivity between viable habitat patches, leading to a rapid loss of biodiversity. We included the number of zero-flow weeks because increases in the frequency of dry-spells directly correlates with reduction in diversity and biomass due to the loss of suitable aquatic habitat. We considered maximum flow because increases in high flow might reduce abundance of young-of-the-year fish by washing away eggs and displacing juveniles and larvae, impeding them from reaching nursery and shelter habitats,.

14We quantified species-specific thresholds for minimum and maximum weekly flow, maximum number of zero flow weeks and maximum and minimum weekly water temperature based on the present-day distribution of these characteristics within the geographic range of each species, similarly to previous studies,,. To this end, we overlaid the species’ range maps with the weekly flow and water temperature metrics from the output of the hydrological model, calculated for each year and averaged over a 30-years historical period to conform to the standard for climate analyses, (1976–2005, for each GCM employed in the study). We calculated for each 5 arcminutes grid cell the long-term average minimum and maximum weekly flow (Qmin, Qmax, Eqs. (1) and (2)), the long-term average frequency of zero-flow weeks (Qzf, Eq. (3)) and the long-term average minimum and maximum weekly temperature (Tmin, Tmax, Eqs. (4) and (5)), as follows:where Q7 and T7 are the vectors of weekly streamflow and water temperature values for a given year i, respectively; q7 is the streamflow value for the week j; N is the number of years considered (30 in this case) and M is the number of weeks in a year (~52). We then used the spatial distributions of these values within the range of each species to determine species-specific ‘thresholds’ for each of the variables, defined as the 2.5 percentile of the minimum flow and minimum temperature and the 97.5 percentile of the maximum water temperature and zero flow weeks values. We preferred these to using the absolute minimum and maximum values to reduce the influence of uncertainties and outliers in the threshold definition. Only for maximum flow we used the maximum value across the range, because of the highly right-skewed distribution of flow values within the range of the species. An overview of the thresholds’ distribution is available in Supplementary Fig. 8.

15We considered four main future scenarios based on increases of global mean air temperature equal to 1.5, 2.0, 3.2, and 4.5 °C. The global mean temperature increase refers to a 30-years average, in accordance with guidelines for climate analyses, and with pre-industrial reference set at 1850–1900. To obtain estimates of weekly water temperature and flow for each warming level, we forced the hydrological model with the output from an ensemble of five Global Climate Models (GCMs), each run for four Representative Concentration Pathway (RCP) scenarios, namely RCP 2.6, 4.5, 6.0, and 8.5 (see “Supplementary Methods” for details). Hence, each RCP–GCM combination would reach each warming level at a different point in time, with some of the RCP–GCM combinations not reaching certain warming levels. Consequently, the number of scenarios available differed among warming levels (an overview is provided in Supplementary Table 1). In total we modeled 42 scenarios (one scenario = one GCM–RCP combination at a certain point in the future), including 17 scenarios for 1.5 °C, 15 for 2.0 °C, 7 for 3.2 °C and 3 for 4.5 °C.

16For each species and each of the 42 scenarios as described in the previous section, we quantified the proportion of the range where projected extremes exceed the present-day values within the species’ range for at least one of the variables. Thus, for each species x we quantified the percentage of geographic range threatened (RT [%]) at each GCM-RCP scenario combination c and for a variable (or group of variables) v as,where AT is the portion of area threatened [km2] and A is the current geographic range size [km2]. That is, we assessed for all grid cells within the species’ range if a projected minimum or maximum weekly flow would fall below the minimum or above the maximum flow threshold, if there would be a higher number of zero flow weeks than the threshold would allow, or if the minimum or maximum weekly water temperature would be lower than the minimum or higher than the maximum water temperature threshold. The variable-by-variable evaluation allowed us to identify which (groups of) variable(s) contributed to the threat. For simplicity, we grouped the number of zero flow weeks, minimum and maximum weekly flow variables to assess threat imposed by altered flow regimes. Similarly, we grouped threats imposed by amplified minimum and maximum weekly water temperature to assess temperature-related threats. In the aggregated results, a grid-cell is thus flagged as threatened if any of the underlying thresholds is exceeded.

17In general, organisms may adapt to climate change (or escape from future extremes) by moving to more suitable locations. Accounting for this possibility is challenging due to the uncertainties and data gaps associated with current and future barriers in freshwater systems (e.g., dams, weirs, culverts, sluices). In addition, data needed to reliably estimate dispersal ability is still lacking for the majority of the species. We therefore employed two relatively simple dispersal assumptions in our calculations. Under the “no dispersal” assumption, fishes are restricted to their current geographic range, whereas under the “maximal dispersal” assumption, fishes are assumed to be able to reach any cell within the sub-basin units encompassing their current geograhic range. We defined the sub-basin units by intersecting the physical boundaries of main basins (defined as having an outlet to the sea/internal sink) with the boundaries defined by the freshwater ecoregions of the world, which provide intra-basins divisions based on evolutionary history and additional ecological factors relevant to freshwater fishes (Supplementary Fig. 11). Basins smaller than 1,000 km2 were combined with adjacent larger units. In total, we delineated 6,525 sub-basin units (area: µ = 20,376 km2, σ = 90,717 km2) from 10,884 main hydrologic basins and 449 freshwater ecoregions. To model future climate threats under the maximal dispersal assumption, we first expanded the geographic range for the current situation, allowing the species to occupy grid cells within the encompassing sub-basin boundaries if suitable according to the species-specific thresholds. Then we assessed future climate threats for the 42 different scenarios relative to the present-day range plus all cells potentially available to the species within the encompassing sub-basins (excluding cells that would become threatened in the future), aswhere AE is the expanded part of the geographic range [km2] and AET is the area threatened within the expanded part of the geographic range [km2].

18To summarize our results, we first assessed the proportion of species having more than half of their (expanded) geographic range threatened (i.e., exposed to climate extremes beyond current levels within their range) at each warming level. We did this for each GCM-RCP scenario combination and then calculated the mean and standard deviation across the GCM-RCP combinations at each warming level. We further calculated the proportion of species threatened by future climate extremes in each 5 arcminutes (~10 km at the Equator) grid cell for each warming level, as follows:where PAF represents the potentially affected fraction of species in grid cell i for warming level w, c represents the scenario (i.e., GCM–RCP combination), Si,w represents the number of species for which extremes in water temperature and flow in grid-cell i according to warming level w do not exceed present-day levels within their range, and Si,present represents the number of species in grid cell i. For both numerator and denominator, the species pool for cell i was determined based on the overlap with the (expanded) geographic range maps (see “Species occurrence data” and “Accounting for dispersal”). We used the median across the GCM–RCP combinations rather than the mean because the data showed skewed distributions. Finally, we averaged the grid-specific proportions of species affected over main basins with an outlet to the ocean/sea or internal sink (e.g., lake), as follows:where Ix represents the number of grid cells within the watershed x.

19We performed phylogenetic regression to relate the threat level of each species, quantified as the proportion of the geographic range exposed to future climate extremes beyond current levels within the range (see Eqs. (6) and (7)), to a number of potentially relevant species characteristics, while accounting for the non-independence of observations due to phylogenetic relatedness among species. We established a phylogenetic regression model per warming level and dispersal scenario (i.e., eight models in total, based on four warming levels times two dispersal assumptions). As species characteristics, we included initial range size (in km2), body length (in cm), climate zone, trophic group and habitat type, as these traits may influence species’ responses to (anthropogenic) environmental change,,,. We further included IUCN Red List category to evaluate the extent to which current threat status is indicative of potential impacts of future climate change, and commercial importance to evaluate implications of potential extirpations for fisheries. We overlaid each species’ geographic range with the historic Köppen–Geiger climate categories to obtain the main climate zone per species (i.e., capital letter of the climate classification). Species falling into multiple climate categories were assigned the climate zone with the largest overlap. We retrieved information on threat status from IUCN and on taxonomy from Fishbase. We used the IUCN and Fishbase data also to gather a list of potential habitats for each species. For the species represented within the IUCN dataset, we classified species as lotic if they were associated with habitats containing at least one of the words “river”, “stream”, “creek”, “canal”, “channel”, “delta”, “estuaries”, and as lentic if the habitat descriptions contained at least one of the words “lake”,“pool”,“bog”,“swamp”,“pond”. For the remaining species, we extracted information on habitat from Fishbase, where we used the highest level of aggregation of habitat types to classify species found in lakes as lentic and species found in rivers as lotic. We classified species occurring in both streams and lakes as lotic-lentic and labeled species found in both freshwater and marine environments as lotic-marine. Further, we retrieved data from Fishbase on maximum body length and commercial importance. From the same database we also retrieved trophic level values and aggregated them into Carnivore (trophic level >2.79), Omnivore (2.19 < trophic level ≤ 2.79) and Herbivore (trophic level ≤ 2.19). We performed a synonym check for the binomial nomenclature provided by the IUCN database to maximize the overlap with the Fishbase database. Since information on phylogeny was available only for a subset of 4,930 fish species covered in our study (based on Betancur-R et al.), we allocated the remaining species to the phylogenetic tree using an imputation procedure implemented in the R package “fishtree”. The empirical tree covered 97% of the families and 80% of the genera included in our species set, indicating that the majority of the missing species were allocated to the correct genus. Our final dataset for the regression included 9,779 species (695 species were excluded because covariates were not available and 951 because they were not included in the “fishtree” database). To account for the uncertainty in the phylogenetic tree imputation, we repeated each of the eight phylogenetic regression models based on 100 replicates of the phylogenetic tree. Prior to running the regressions we log-transformed threat level (response variable), geographic range size and body length as these variables were right-skewed. As Spearman’s rank correlations among the covariates were below 0.4 and variance inflation factors below 1.5 (Supplementary Fig. 12 and Supplementary Table 6), we kept the full set of covariates. We ran the phylogenetic regression using the R package “nlme”, and checked the residuals of the models using QQ plots (Supplementary Fig. 13). We then extracted coefficients, |t|-statistics, p values as well as the lambda parameter at each warming level (averaged over the 100 replicates; Supplementary Table 5). Then, we quantified variable importance using a procedure based on the random forest approach, as implemented in the R package “biomod2”. To that end, we randomized the values of the covariates one by one and computed the variable importance as 1 minus the Pearson’s r between the predictions of the original model and the predictions obtained from the model with randomized data. We iterated this procedure 10,000 times (100 iterations of the variable importance algorithm times 100 models based on replicates of the phylogenetic tree) and reported the average score and standard deviation across the 100 stochastic replicates (standard deviation across the iterations was negligible) for each of the eight models.

20Further information on research design is available in the Nature Research Reporting Summary linked to this article.