Linking environmental salinity to respiratory phenotypes and metabolic rate in fishes: a data mining and modelling approach

Seawater (SW) and freshwater (FW) environments present different osmoregulatory challenges to resident species. Teleosts in the two environments have similar internal ion concentrations (Holmes and Donaldson, 1969) and, therefore, are hypo-osmotic compared with ion-rich SW, where they risk losing water to the environment, and hyper-osmotic compared with ion-poor FW, where they risk losing ions (Evans et al., 2005). The fish gill is the primary site of gas exchange and presents a large, thin epithelium to the surrounding water that facilitates the uptake of oxygen (O₂) and the excretion of the metabolic by-products carbon dioxide and ammonia. However, those gill morphometrics that maximize diffusive gas exchange may also increase passive ion and water fluxes that, if unopposed, would lead to osmotic disturbances. These conflicting requirements for ionoregulation and O₂ uptake at the fish gill were recognized over 100 years ago by August Krogh (Ege and Krogh, 1914; Krogh, 1937) and later substantiated in seminal work by Randall et al. (1972), Gonzalez and McDonald (1992) and others. Today, the effect is known as the osmorespiratory compromise, and it has been treated in detailed reviews (Gilmour and Perry, 2018; Gonzalez, 2011; Nilsson et al., 2012; Perry, 1998; Sardella and Brauner, 2007; Wood and Eom, 2021). The present Review builds on previous work to uncover general patterns among teleost fishes that link environmental salinity to gill and blood phenotypes and to provide mechanistic insight into the physiology underlying these observations by using a modelling approach.

FW and SW teleosts have evolved divergent strategies to overcome the osmotic challenges imposed by their respective environments, involving anatomical and physiological adaptations at different organ systems. Briefly, FW fishes take up ions at the gills and excrete dilute urine via the kidneys, whereas SW fishes take up water and ions via their digestive tracts and excrete ions via the gills and kidneys, and with their faeces (Evans et al., 2005; Grosell, 2010; Marshall and Grosell, 2006). These divergent mechanisms enable FW and SW fishes to maintain water and ion homeostasis in their respective environments and are well studied, but their effects on gill morphometrics and, thus, O₂ uptake capacity, have not been studied broadly. To address this knowledge gap, we conducted a systematic literature review of gill, blood and metabolic characteristics in 355 teleost species in relation to their aquatic habitat (155 FW, 200 SW species; see Appendix for details on the literature review and Table S2 for raw data).

Our initial analysis indicated that FW fishes, on average, had significantly smaller mass-specific gill surface area (GSA) than SW fishes (P=0.019; Fig. 1A), consistent with what has long been reported in the literature. However, we also observed differences in body mass between the investigated FW and SW fishes, which biases this analysis because GSA scales allometrically with body mass. Furthermore, these interspecific comparisons are dominated by a few hyper-diverse taxa (Ostariophysi in FW and Percomorpha in SW; Vega and Wiens, 2012), which may introduce an additional phylogenetic bias. Therefore, to test whether salinity affected GSA while taking into account allometric scaling and phylogenetic non-independence of species, we fitted a phylogenetically corrected linear model to the log₁₀-transformed GSA and body mass data (Fig. 1B), and then calculated the vertical residuals that provided a body mass-independent measure of GSA (Garland et al., 1992). These residuals were independent of salinity (P=0.672; Fig. 1C), indicating that, once corrected for mass and phylogenetic relatedness, teleosts in FW and SW had similar GSAs. These findings illustrate the importance of accounting for body mass and species relatedness in interspecific comparisons of traits that scale allometrically with body mass. In contrast, the blood–water diffusion distance (hereafter, gill thickness; Fig. 1D) was unaffected by body mass (P=0.789; Fig. 1E) and, thus, these data were only corrected for the phylogenetic relatedness of the species. As with GSA, we found no significant effect of salinity on gill thickness (P=0.140; Fig. 1F), indicating that other factors must determine the observed variability.

To explore the interacting effects of GSA and thickness, we calculated the O₂ diffusive conductance of the gills (G_d) in FW and SW species, which describes the amount of O₂ that can diffuse across the gills for a given P_O2 gradient. In our study, G_d was calculated from GSA and gill thickness values from the literature, and from the gill characteristics described in Table S1 (Dejours, 1981). We found 33 species with reported values for both GSA and thickness (Fig. 1G) and, after correcting for body mass and phylogenetic relatedness (Fig. 1H), FW fishes had significantly lower G_d values than SW fishes (P=0.017; Fig. 1I). This result should be interpreted in light of the relatively few species with thickness measurements and the fact that some species with extreme gill characteristics may disproportionately influence the analysis (such as the highly active Scombridae in SW and some facultative air-breathing fishes in FW). We tried expanding our analysis to all actinopterygian fishes to obtain additional values for G_d but found that gill thickness values were also largely unavailable. Additional measurements of gill thickness in other FW and SW fishes are needed to strengthen the relationships between G_d and environmental salinity that we found in our analysis. Nevertheless, the data that are currently available indicate that, even after mass and phylogenetic corrections, FW teleosts have lower G_d than SW teleosts, and if substantiated more broadly, this finding may change our understanding of the physiological implications of environmental salinity on gill phenotypes.

Smaller and thicker gills can benefit fishes by reducing the rate of passive ion and water fluxes that can disrupt osmotic homeostasis (Greco et al., 1996; Henriksson et al., 2008; Perry, 1998; Sollid et al., 2003). However, the gill is also the primary site for gas exchange and smaller and thicker gills will result in lower G_d. In FW fishes, lower G_d may limit O₂-uptake capacity and could impair their performance in situations of high O₂ demand (e.g. exercise) or low O₂ availability (e.g. environmental hypoxia). Interestingly, our data revealed a potential compensatory mechanism, as FW species had a ∼2-fold higher haemoglobin (Hb) O₂ affinity than SW species (expressed as lower P₅₀ values of Hb – the P_O2 at which Hb is 50% saturated with O₂; P<0.001; Fig. 1L). Hb P₅₀ values were corrected for phylogenetic relatedness and a significant temperature effect using a Van't Hoff plot (P = 0.001; Fig. 1K). Because Hb is the principal O₂ carrier in the blood and its O₂-binding characteristics determine the kinetics of O₂ loading at the gills and unloading at the tissues, the higher Hb–O₂ affinity of FW species could enhance the diffusion gradient of O₂ across the fish gill and theoretically support higher rates of metabolic O₂ consumption (Ṁ_O2; a proxy for metabolic rate). These observations led us to hypothesize that the high Hb–O₂ affinities of FW species compensate for their lower G_d, allowing FW species to maintain similar Ṁ_O2 as SW species during aerobic challenges, such as exercise and hypoxia.

To test this hypothesis, we explored the interacting effects of G_d and Hb P₅₀ using a mathematical model (based on Malte and Weber, 1985) that predicts the maximal Ṁ_O2 (Ṁ_O2,max) that a fish can support over a range of gill and blood characteristics, as well as environmental conditions (see Appendix for detailed model description). Specifically, we tested three predictions derived from our hypothesis: all else being equal, (1) a higher Hb–O₂ affinity will compensate for a lower G_d to maintain arterial O₂ transport; (2) FW and SW teleosts will achieve similar Ṁ_O2,max using different combinations of G_d and Hb–O₂ affinity; and (3) FW and SW teleosts will achieve similar hypoxia tolerance (represented by P_crit – the lowest water P_O2 at which a fish can maintain its standard Ṁ_O2; Ṁ_O2,std) using different combinations of G_d and Hb–O₂ affinity. We then compared the modelled results with empirical values mined from the literature that revealed novel patterns in the divergent gill and Hb phenotypes of FW and SW teleosts in relation to environmental salinity.

First, we used our Ṁ_O2,max model to investigate the effects of gill and blood characteristics on the equilibration of P_O2 across the gill epithelium using the lowest, median and highest G_d values for FW and SW fishes that we observed in the literature (using only values that were based on N>1), over the range of Hb–O₂ affinities reported in Fig. 1. FW fishes with the lowest G_d were severely diffusion-limited, reflected in a poor equilibration of P_O2 between the counter-current flows of water and blood (Fig. 2A). Ṁ_O2,max in these fishes was generally low, but the highest values were achieved at the lowest Hb P₅₀. Therefore, when all else is equal, a higher Hb–O₂ affinity can partially compensate for a low G_d by increasing the P_O2 gradient across the epithelium. This finding is consistent with previous work that modelled the diffusion of O₂ across the fish gill and found an increased O₂ extraction from the water and lower ventilatory requirements at high Hb–O₂ affinities (Malte and Weber, 1987). In addition, a reduction in Hb–O₂ affinity caused a steep decline in the arterial O₂ saturation of Hb (Sa_O2), which determines the maximum capacity for blood O₂ transport and Ṁ_O2,max (Gallaugher et al., 2001). In SW species, the lowest, median and highest values for G_d were higher than in FW fishes. Consequently, O₂ uptake at the gills in SW was generally more efficient, as indicated by a nearly complete equilibration of P_O2 across the gill epithelium, especially at the lowest Hb P₅₀ values (Fig. 2D).

An increase in G_d to median values increased Ṁ_O2,max in both FW and SW fishes and shifted the optimal Hb P₅₀ to higher values (Fig. 2B,E). Higher G_d also reduced the effect of Hb P₅₀ on Sa_O2 because O₂ uptake became less dependent on the P_O2 diffusion gradient. The fact that a high Hb–O₂ affinity exerts the greatest benefits in fish with the lowest G_d is consistent with our literature-mined data. Combined, these findings indicate that the higher Hb–O₂ affinities observed in FW fishes may have adaptive significance to compensate for their low-G_d gills, an idea that could be substantiated by investigating G_d in additional species.

SW fishes with median G_d achieved full O₂ saturation during blood transit through the gill, suggesting that Ṁ_O2,max is limited by maximal cardiac output (⁠⁠), rather than the diffusion characteristics of the gill (Fig. 2D), and the same was true for those FW and SW species with the highest G_d (Fig. 2C,F). These findings generally agree with experimental work, indicating that O₂ uptake at the gill is not diffusion-limited in normoxia, but rather is a function of maximal perfusion and/or ventilation rates (Daxboeck et al., 1982; Malte and Weber, 1985; Randall and Daxboeck, 1984). However, in >50% of the FW species in our mined dataset, Ṁ_O2,max may be limited by gill diffusion characteristics, where increases in are inconsequential. These results must be interpreted with respect to the value for rainbow trout that we used in our simulations, which may be higher than those of many low-G_d FW species. Nevertheless, our results highlight systematic differences in cardiorespiratory phenotypes between FW and SW teleosts that have ecological consequences. Those species with the highest G_d may achieve increases in Ṁ_O2,max through increases in ⁠, enabling more active lifestyles that are not available to low-G_d species. A broader assessment of G_d, Ṁ_O2,max and in relation to environmental salinity in teleosts may shed light on the mechanistic basis for our observations and their ecological significance.

Finally, our simulations also showed that the Hb P₅₀ values that achieved the highest Sa_O2 did not necessarily produce the highest Ṁ_O2,max. Lower Hb P₅₀ values can safeguard Sa_O2 and improve Ṁ_O2,max in those fish with the lowest G_d. However, the higher G_d we observed in SW species would lessen the benefits of a lower Hb P₅₀, potentially enabling them to exploit a higher range of Hb P₅₀ values that promote the offloading of O₂ at the tissues, which we explore in more detail below.

The optimal Hb P₅₀ is a compromise between the physiological requirements for O₂ loading at the gas exchange surface and unloading at the tissue capillaries (Brauner and Wang, 1997; Harter and Brauner, 2017; Wang and Malte, 2011). To highlight this trade-off in more detail, Fig. 3 shows the Ṁ_O2,max that can theoretically be attained over the range of Hb P₅₀ and G_d values in FW and SW teleosts. Generally, those species with the lowest G_d achieved the lowest Ṁ_O2,max and increasing G_d improved Ṁ_O2,max, consistent with a diffusion limitation on O₂ uptake. However, Ṁ_O2,max in species with median G_d matched those of species with the highest G_d, indicating a progressive perfusion limitation, which is in line with the results shown in Fig. 2. As the upper boundary for Ṁ_O2,max was set by perfusion, the two groups achieved similar Ṁ_O2,max of ∼13 µmol g⁻¹ h⁻¹ at the for rainbow trout (53 ml kg⁻¹ min⁻¹).

Increasing Hb–O₂ affinity also had beneficial effects on Ṁ_O2,max, but only until an optimum value was reached, after which further increases in Hb–O₂ affinity severely decreased Ṁ_O2,max. The effect of Hb–O₂ affinity on Ṁ_O2,max depended on G_d, where a reduced G_d shifted the optimal P₅₀ for Ṁ_O2,max to lower values. For example, when G_d was reduced from median to lowest values, the optimal P₅₀ shifted from 31 to 3 mmHg in FW fishes (Fig. 3A) and from 39 to 25 mmHg in SW fishes (Fig. 3B). These ranges of optimal P₅₀ map well onto the mined P₅₀ values for FW and SW fishes, respectively (Fig. 1).

Our simulations are also in line with data on rainbow trout acclimated to soft water that thickened their gills from 3 to 6 µm as a result of ionocyte proliferation and compensated for the impaired G_d by decreasing P₅₀ from 18 to 12 mmHg (Greco et al., 1996; Perry, 1998; Perry et al., 1996). Gill morphology and P₅₀ are plastic traits, and it seems that they may respond in concert to overcome physiological challenges to O₂ uptake in a way that depends on environmental salinity. The situation is complicated as the diffusive pathways for O₂, water and ions may differ and many species exert some control over these ‘passive’ fluxes, such as observed in extremely hypoxia-tolerant species (Wood et al., 2009), active species (Gonzalez and Mcdonald, 1994), and after exercise training (Gallaugher et al., 2001; Postlethwaite and McDonald, 1995). A salinity-specific effect on gill permeability has not been studied broadly, but data for a few species exist. FW-acclimated killifish have higher branchial water permeability than SW conspecifics, but they display similar diffusive ion fluxes (in opposite directions) during hypoxic hyperventilation (Giacomin et al., 2019; Wood et al., 2019). However, while FW-acclimated killifish do not alter ionocyte density, they do decrease their overall GSA (Giacomin et al., 2019), which is consistent with our observation of lower G_d in FW fishes that may perhaps occur via different mechanisms in different species and conditions. If this effect is substantiated more broadly, one may hypothesize that the hyper-osmoregulatory strategy of FW fishes has led not only to lower G_d but also to higher Hb–O₂ affinities that balance the requirements for O₂ uptake and tissue O₂ extraction at the prevailing gill diffusion characteristics. Likewise, the hypo-osmoregulatory strategy of SW fishes may be permissive of high-G_d gills that are best matched by higher Hb P₅₀ and that enable increases in and Ṁ_O2,max in species with active lifestyles where exercise performance is linked to evolutionary fitness.

To further explore the interacting effects of G_d and Hb–O₂ affinity on Ṁ_O2,max, we plotted the upper and lower boundaries for G_d as a function of Hb P₅₀, resulting in areas of theoretical Ṁ_O2,max that can be achieved by the gill and blood characteristics of FW and SW fishes (Fig. 4A). The lower areas were calculated at the for rainbow trout, and because of the perfusion limitation of Ṁ_O2,max, FW and SW fishes with the highest G_d achieved similar Ṁ_O2,max. To corroborate these modelling results, we then plotted empirically measured Ṁ_O2,max values onto the theoretical areas. Generally, the measured and calculated values agreed well. However, one cluster of FW fishes with high Hb–O₂ affinity fell above the predicted lines, which is likely explained by higher tissue O₂ conductance that maintains tissue O₂ extraction from high-affinity Hbs, and perhaps by lower Hill coefficients that enhance O₂ extraction from the water (Malte and Weber, 1987); both these values were kept constant in our model. Other FW and SW species may achieve values that exceed those of rainbow trout. For instance, sockeye salmon swimming in FW reached Ṁ_O2,max of 27.2 µmol g⁻¹ h⁻¹ at a of 67.8 ml kg⁻¹ min⁻¹ (Steinhausen et al., 2008) and in some sockeye salmon populations with longer migration routes, can exceed 100 ml kg⁻¹ min⁻¹, nearly 2-fold the value of rainbow trout (Eliason et al., 2013). Whether the gill characteristics of these anadromous sockeye salmon are truly representative of a FW species is unclear, but despite their semelparous life history, they appear to exert tight control over ion homeostasis, even during their final FW stage.

We then explored whether a salinity-specific pattern in Ṁ_O2,max exists across empirically measured values for 130 fish species (51 FW, 79 SW), and we corrected these data for allometric scaling with body mass, temperature and phylogenetic relatedness (Fig. 4B). Across all species, there was no significant effect of salinity on Ṁ_O2,max (P<0.067; Fig. 4C), consistent with our prediction that FW and SW teleosts may achieve similar Ṁ_O2,max through different combinations of gill and blood characteristics. However, the low P-value in our analysis on 130 species suggests that future studies, especially those with larger sample sizes, may be able to resolve potential differences in Ṁ_O2,max between FW and SW teleosts. Finally, Ṁ_O2,max is a complex physiological trait that is strongly influenced by methodology, animal condition and behaviour (Killen et al., 2017; Norin and Clark, 2016); exploring these sources of variability and accounting for them by statistical means may be a worthwhile avenue for future analyses.

Some of the highest Ṁ_O2,max values are found in the fast-swimming, pelagic, SW teleosts (Wegner et al., 2010), including billfishes (Istiophoridae and Xiphiidae), dolphinfishes (Coryphaenidae), jacks (Carangidae) and tunas (Scombridae). Reliable data during maximal exercise in these animals are notoriously difficult to obtain, but a few measurements are available. Based on submaximal swimming trials, the Ṁ_O2,max in skipjack tuna has been estimated at 68.7 µmol g⁻¹ h⁻¹ (Brill, 1987; Dewar and Graham, 1994; Gooding et al., 1981) and therefore exceeds the value for rainbow trout by more than 5-fold. These high Ṁ_O2,max values are enabled by GSAs of nearly 20 cm² g⁻¹ in yellowfin and skipjack tunas (Brill and Bushnell, 2001), and gill thicknesses as low as 0.5 µm (Hughes, 1970, 1984), resulting in G_d values that are 3-fold higher than the highest values of any FW species. Clearly, a unique physiology in tunas is required to sustain such high Ṁ_O2,max that cannot be achieved with the cardio-respiratory characteristics of rainbow trout or most other teleosts. Therefore, to determine what physiological adaptations are required in tuna to attain their high Ṁ_O2,max, we re-calibrated our model to the skipjack tuna (Fig. 4A; light-shaded areas). This adjustment revealed that the high Ṁ_O2,max in tuna is feasible only if is increased ∼6-fold (to 318 ml kg⁻¹ min⁻¹), blood [Hb] is increased 2.3-fold and O₂ conductance at the tissues is increased ∼3-fold over the values in rainbow trout. Whether these values are representative of some tuna species remains to be validated experimentally; however, the available data indicate that they may not be unreasonable estimates. Even in spinally blocked skipjack tuna, has been measured at 132 ml kg⁻¹ min⁻¹ (Brill, 1987; Bushnell, 1988), blood [Hb] can reach 2.3 mmol l⁻¹ (Brill and Bushnell, 1991) and tissue O₂ conductance is undoubtedly increased by extremely high capillary density (Hulbert et al., 1979) and myoglobin concentrations (George and Stevens, 1978; Stevens and Carey, 1981). Theoretically, high ⁠, [Hb] and tissue O₂ conductance may also be attainable by FW teleosts; but matching high-G_d gills may not. Our simulations show that even the highest G_d values in FW fishes are insufficient to produce tuna-like Ṁ_O2,max and the upper boundary levels off at ∼40 µmol g⁻¹ h⁻¹. These values are still higher than any Ṁ_O2,max recorded for a FW teleost, indicating that factors other than G_d may set the actual limit. Regardless, these fundamental relationships may explain why no FW fishes match the high Ṁ_O2,max found in some SW species. Whether FW tuna-like fishes are, in fact, absent from the fossil record is a worthwhile avenue for future investigations. However, based on our findings it would seem that FW environments should remain devoid of fast-swimming pelagic analogues to the tunas.

Another important driver for variation in gill characteristics and Hb–O₂ affinity among fishes is the availability of environmental O₂ (Mandic et al., 2009). Under conditions of low environmental O₂ (hypoxia), a fish's indefinite survival depends on its ability to sustain Ṁ_O2,std (Hughes, 1973). The lowest water P_O2 at which a fish can sustain Ṁ_O2,std is termed the critical O₂ tension (P_crit), and at P_O2 below P_crit, a fish's survival becomes dependent on some combination of anaerobic glycolysis and metabolic depression, both of which are unsustainable in the long term (Beamish, 1964; Regan et al., 2017; Ultsch and Regan, 2019; Ultsch et al., 1978). P_crit is therefore a common metric of hypoxia tolerance that represents the suite of aerobic contributions to hypoxia survival in a single value (Regan et al., 2019); thus, we explored this metric in our simulations.

To estimate how differences in G_d and Hb–O₂ affinity may influence P_crit, we adjusted our model to predict Ṁ_O2,max at progressively decreasing inspired water P_O2 (from 150 to 0 mmHg) and then fixed a horizontal line representing Ṁ_O2,std at 1.5 µmol g⁻¹ h⁻¹, which is the value reported for resting rainbow trout (Kiceniuk and Jones, 1977). The P_O2 at which the Ṁ_O2,max and Ṁ_O2,std lines intersect represents the theoretical P_crit for a given set of gill and blood characteristics. Our approach assumed that Ṁ_O2,max and Ṁ_O2,std respond similarly to reductions in P_O2 below P_crit, and that the typical P_crit curve is biphasic. The former assumption is consistent with previous theoretical work (Esbaugh et al., 2021), empirical measurements (Claireaux et al., 2000) and the fact that aerobic scope for activity at sub-P_critP_O2 is zero (Claireaux and Chabot, 2016; Ern et al., 2016). The latter assumption is true for some fish species, but not all (Wood, 2018; Farrell et al., 2021). However, most important for our simulations was not the shape of the Ṁ_O2,std curve at P_O2 below P_crit, but rather the presence of a benchmark Ṁ_O2,std value to which the sub-P_critṀ_O2,max line could be compared. A biphasic P_crit curve enabled this.

The outcomes of these simulations are shown in Fig. 5 for the lowest, median and highest G_d that we observed in FW and SW fishes (Fig. 1). In all cases, increasing Hb–O₂ affinity or G_d led to lower P_crit values, representing a higher hypoxia tolerance of the fish. SW fishes achieved lower P_crit than FW fishes when G_d was at the lowest or median values. However, the lowest overall P_crit values were achieved by FW fishes with high G_d and high Hb–O₂ affinity. These and our previous results show that an increased Hb–O₂ affinity can generally compensate for impaired diffusion characteristics at the gill caused by either low G_d or low environmental O₂. In both cases, Ṁ_O2 is limited by the reduction in Sa_O2 that can be improved by increasing Hb–O₂ affinity, and these results are in close agreement with previous work (Wang and Malte, 2011).

To validate our simulations, we plotted the P_crit values for high and low G_d as a function of Hb–O₂ affinity, spanning the area of theoretical P_crit that were predicted by the model. We then overlayed empirically measured P_crit values from the literature for FW and SW fishes (Fig. 6A). For FW fishes, we found an excellent agreement between the model predictions and measured P_crit values. In addition, P_crit values in FW fishes were highly variable, spanning an order of magnitude. GSA and gill thickness are plastic traits that respond to hypoxia, environmental ion concentration and temperature. The underlying mechanisms include dynamic changes in the redistribution of branchial blood flow and cellular responses that reversibly alter the morphology of the gill (Wood and Eom, 2021). For example, the crucian carp can increase GSA by ∼7.5-fold after an acclimation period to hypoxia, largely through the receding of an interlamellar cell mass (Sollid et al., 2003). These acclimation responses can occur quickly, as in goldfish, where GSA increases ∼2-fold in just hours (Regan et al., 2017). High GSA plasticity is often found in hypoxia-tolerant rather than -intolerant species (Dhillon et al., 2013), and in more FW than SW species (Gilmour and Perry, 2018). This may be because many FW environments are particularly hypoxia-prone (Diaz and Breitburg, 2009) or simply because relatively few SW species have been studied in this respect.

In contrast, SW fishes occupied a narrow band in P_crit values over a broad range of Hb P₅₀, which is supported by the modelled and empirical data (Fig. 6A). However, the modelled data were generally offset to lower P_crit values. The reasons for this discrepancy may relate to our model assumptions not accurately representing some of the in vivo characteristics of SW fishes; for example: (i) many SW species may have higher Ṁ_O2,std than rainbow trout, and/or (ii) the boundaries for G_d in our dataset may not match those of SW species for which P_crit and Hb P₅₀ data are available, both representing only a small subset of SW species. Regardless, the combined data revealed that FW and SW fishes occupy different quadrants in the P_crit versus Hb P₅₀ relationship, which indicates that they achieve their hypoxia tolerance by different mechanisms. SW fishes rely on their higher G_d to maintain branchial O₂ uptake in hypoxia (Fig. 5), resulting in P_crit values that are largely independent of Hb P₅₀ and thus fall within a narrow range (Fig. 6A). In contrast, FW fishes rely on their high Hb–O₂ affinity to maintain branchial O₂ uptake in hypoxia and, when combined with high G_d values, achieve the lowest P_crit observed in teleosts.

To test whether there is a systematic difference in P_crit between FW and SW teleosts, we mined literature values for 142 species (46 FW and 96 SW). After correcting for phylogenetic relatedness of species, we found no significant difference in P_crit between FW and SW teleosts (P=0.891; Fig. 6C). This stands in contrast with our dataset when not corrected for phylogeny, which showed significantly lower P_crit in FW than in SW fish, and to the metanalysis of Rogers et al. (2016) and intraspecific salinity trials of Haney and Nordlie (1997), which reported lower P_crit in FW fishes under some, but not all, conditions. Therefore, the variation in P_crit across salinity environments appears to be driven by a few teleost clades with extreme P_crit values, again highlighting the importance of applying appropriate corrections of the raw data. Relatedly, most SW species on which P_crit experiments have been performed are native to environments that regularly experience hypoxia, such as coral reefs and intertidal zones. The mean SW P_crit value may therefore be skewed towards these hypoxia-adapted SW species and away from pelagic species that are less likely to regularly encounter hypoxia. Fewer P_crit values for pelagic species exist, perhaps because these fishes are not easily caught and experimented on. However, their addition to this comparison could further distinguish the salinity groups and is a worthwhile direction for future studies. Finally, there is substantial variability in P_crit values in the literature and, as with other complex performance metrics, P_crit should be interpreted with some caution regarding differences in methodology, animal condition and behaviour (Wood, 2018). Part of this variability may be due to the plasticity in gill morphology and blood characteristics in teleosts, perhaps coupled with different methodologies regarding the rates of hypoxia onset (Rogers et al., 2016). Nevertheless, the close agreement of the empirical and modelling data indicates that the general patterns in P_crit between FW and SW fishes are largely driven by factors that are manipulated in the model. Therefore, studying the interspecific variability in G_d, Hb–O₂ affinity and their interactions may provide a roadmap for future investigations into the mechanistic basis for P_crit in fish.

Our analysis provides compelling evidence for systematic differences in gill and blood phenotypes between FW and SW teleosts that may be related to the osmotic characteristics of their aquatic environments. However, aerobic performance, assessed as Ṁ_O2,max and P_crit, did not differ between FW and SW fishes; therefore, the two groups may use different combinations of gill and blood characteristics to support similar aerobic capacities. FW fishes generally had lower G_d, which has consequences for aerobic metabolism and may reflect a physiological constraint by their hyper-osmoregulatory strategy. The higher Hb–O₂ affinity of FW fishes may thus have adaptive significance, as it increases the diffusion gradient for O₂ across the gills. The benefits of a higher Hb–O₂ affinity may not fully compensate for the disadvantages of a low G_d, but rather, will enable a higher Ṁ_O2,max within this constraint, by balancing the conflicting requirements for O₂ loading at the gills and unloading at the tissues. In contrast, the hypo-osmoregulatory strategy of SW teleosts is clearly permissive of higher G_d, which in turn may have lifted the brakes on the evolution of higher and [Hb]. In tuna, these coordinated adaptations of their cardio-respiratory systems ultimately enabled extremely high Ṁ_O2,max that are not attainable even by FW fishes with the highest G_d.

Gill morphology and Hb–O₂ binding characteristics are highly plastic traits in teleosts, and many species with euryhaline, anadromous or catadromous life cycles routinely transition between the FW and SW environments and similar transitions occurred repeatedly in the course of teleost evolution (Betancur-R et al., 2015). The evolutionary dynamics that govern these transitions are worthy of further investigation and fish clades that include species that have recently transitioned between FW and SW may be well suited for this purpose.

The osmorespiratory compromise describes the conflicting requirements for O₂ uptake and ion/water flux at the gills and has typically been studied with respect to GSA. Our data indicate that gill thickness and G_d are important factors that determine Ṁ_O2,max in fishes and should be considered in future studies on the osmorespiratory compromise by measuring O₂ uptake and unidirectional ion/water flux in more species in both FW and SW. Ideally, future studies should use standardized protocols for determining gill morphometrics, unidirectional ion/water flux, Ṁ_O2,max, P_crit and Hb characteristics in animals acclimated to a clearly defined set of environmental conditions, allowing for more reproducible interspecific comparisons. These future studies may lay the foundation to test our hypothesis that invasions of the FW environment and the divergent selective pressures of the osmorespiratory compromise gave rise to adaptations that reduced G_d and Hb P₅₀, thus profoundly shaping the aerobic pathway of FW fishes.

We conducted literature reviews of five cardio-respiratory variables in teleost fishes: gill surface area (GSA), blood–water diffusion distance (gill thickness), Hb–O₂ affinity (expressed as P₅₀, the P_O2 at which Hb is 50% saturated with O₂), maximal O₂ consumption rate (Ṁ_O2,max) and critical O₂ tension (P_crit, the lowest P_O2 at which a fish can maintain standard Ṁ_O2). We included only teleost species in our datasets, which are by far the most diverse group of fishes and account for ∼95% of all extant species (Nelson et al., 2006). The reasons for excluding other groups, such as the Chondrichthyes, were based on their osmo-conforming physiology that changes the relationship between respiratory and ionoregulatory demands at the gill, and this group is almost exclusively found in the marine environment. We also excluded amphibious and obligatory air-breathing fishes (Damsgaard et al., 2020), as these traits fundamentally change the relationship between osmoregulation and respiration that we sought to examine (Graham, 1997). Additionally, the values we used were limited to those from individuals in juvenile or adult developmental stages, measured under control conditions (i.e. untreated and non-acclimated individuals).

Commenting on the methodological approaches used in literature studies was outside the scope of this Review and, generally, we considered data generated with different methods, while being aware that our variables of interest are sensitive to different experimental protocols. However, we did limit our search to Hb P₅₀ measured on whole blood, and excluded those on haemolysates, which are less representative of in vivo conditions even when standardizing for allosteric effectors (Berenbrink, 2006). For GSA, we limited our search to measurements of total GSA or lamellar surface area.

All literature searches were performed on Clarivate's Web of Science^®, in 2020 and 2021, using the following search terms: GSA (gill AND surface area AND fish; lamellar surface area AND fish); Hb P₅₀ (oxygen* binding affinity AND fish; hemoglobin* AND P₅₀ AND fish); P_crit (critical oxygen* AND fish; P_crit+fish; P_crit+hypoxia+fish); Ṁ_O2,max (maximum metabolic rate+fish; MO2max+fish). In instances where the mined studies cited other studies that our search did not uncover, we included these data if they were appropriate. Two recent reviews were particularly valuable in our search, providing many of the mined data points: Killen et al. (2017) and Rogers et al. (2016). Values, references and URLs for all cited studies can be found in Table S2.

We analysed our data in R Studio (v.1.3.1056) using a generalized linear mixed model (GLMM) to identify correlations between continuous variables using the MCMCglmm() function in the MCMCglmm package in R (Hadfield, 2010; Hadfield and Nakagawa, 2010). This method treats individual studies on the same species as a random effect (i.e. taking into account intraspecific variation) and accounted for phylogenetic non-independence of species. Here, we used a maximum clade credibility tree generated from 100 Bayesian posterior probability trees from Rabosky et al. (2018) using the maxCladeCred() function in the phangorn package in R (Schliep, 2011).

To test the effect of salinity on GSA and G_d, while taking into account the hypo-allometric scaling with body mass, we fitted a GLMM model to the log₁₀-transformed data and calculated the residuals. These residuals provide a mass-independent measure of GSA and G_d (Garland et al., 1992). We then tested for an effect of salinity on the residuals with a phylogenetic analysis of variance simulation (Garland et al., 1993), using the phylANOVA() function in the phytools package in R (Revell, 2012). Similar approaches were used for the other variables, as detailed below.

We selected the full model based on the lowest deviance information criterion value (DIC), calculated the residuals, and tested for an effect of salinity on the Ṁ_O2,max residuals, as described above. To test for the effect of salinity on Hb P₅₀ we first fitted a GLMM model to log₁₀(P₅₀) versus absolute temperature^–1 and calculated the residuals that we tested for an effect of salinity. P_crit and gill thickness were body mass independent; thus, we used the raw values in a phylogenetic ANOVA to test for an effect of salinity.

To study the effects of Hb P₅₀ and G_d on gas exchange at the fish gills, we implemented the previously published mathematical model of Malte and Weber (1985) using R v.3.6.2. in RStudio v.1.2.5033. Briefly, the model solves a system of differential equations to calculate the changes in water and blood P_O2 during counter-current gas exchange along the length of a model fish gill. In R, the differential equations were solved as a boundary value problem with venous P_O2 and inspired water P_O2 as the input parameters, using the bvpsolve package (Mazzia et al., 2014), with bvptwp() that uses a mono-implicit Runge–Kutta (MIRK) method with deferred corrections (Cash and Mazzia, 2005; Cash and Wright, 1991) and a continuation method. The output from these gas exchange simulations gives the calculated values for expired water P_O2 and arterial P_O2, which is then used to calculate other arterial blood parameters. Next, the model predicts the outcome of gas exchange at the tissues to calculate venous P_O2 and other venous blood parameters that are then used for the next iteration of gas exchange at the gills; the cycle is repeated until the final calculated venous P_O2 is in equilibrium with the gas exchange systems at the gills and at the tissues (Malte and Weber, 1985, 1987, 1989, 2011). Our model focused only on the exchange of O₂ at the fish gill and did not take into account changes in CO₂ excretion; however, the model does simulate a pH shift at the tissues that leads to an increase in Hb P₅₀ by the Bohr effect and recovery of Hb P₅₀ during venous transit.

The model was ‘calibrated’ against the cardio-respiratory characteristics measured in rainbow trout (Oncorhynchus mykiss) at 10°C by Kiceniuk and Jones (1977). Specifically, the diffusion coefficient for O₂ at the gills (D_O2) and the tissue O₂ conductance, two values for which no reliable measurements exist, were adjusted so that the model output matched the empirically measured Ṁ_O2,max, arterial and venous P_O2, and tissue O₂ extraction in rainbow trout. All other input parameters for the model calibration are listed in Table S1 with their respective references.

In all simulations, [Hb] was set to 1 mmol l⁻¹ based on the summary data by Gallaugher and Farrell (1998) and the cooperativity of Hb–O₂ binding (the Hill coefficient) was set to n=2. While this represents a simplistic approach, the effects of cooperativity on gas exchange at the fish gill have been addressed previously (Malte and Weber, 1987) and we observed no systematic differences in cooperativity between FW and SW teleosts that would warrant a more detailed investigation within the context of our study. In addition, some groups (e.g. Cottidae) display a wide range of Hb P₅₀ at rather constant Hill coefficients, indicating that phylogeny may be a better predictor of Hill coefficients than, for example, Hb P₅₀ (Mandic et al., 2009). Finally, to test our hypotheses, we adjusted the parameters of G_d and Hb P₅₀ based on the values for FW and SW fishes that resulted from our review of the literature.

We thank Nick Wegner and Rod Wilson for insightful discussions and help with our literature review, Tommy Norin for critical feedback on our manuscript, and the two anonymous reviewers for their insightful comments.