Interspecific variation in hypoxia tolerance and hypoxia acclimation responses in killifish from the family Fundulidae

Hypoxia is a pervasive stressor in aquatic environments that can originate from natural and anthropogenic events (Breitburg et al., 2009; Diaz and Breitburg, 2009; Diaz and Rosenberg, 2008). Aquatic hypoxia is most simply defined as a reduction in water O₂ levels that compromises some aspect of physiological function (Farrell and Richards, 2009). Severe bouts of hypoxia are often implicated in fish kills, and hypoxic episodes can have broad ecological implications by reducing available habitat, altering species distributions, and changing trophic relationships (Breitburg et al., 2009; Mallin et al., 2006). Even brief exposure to less severe reductions in O₂ availability can also affect fish physiology. Hypoxia-prone zones like tidal pools and estuaries are typically occupied by relatively hypoxia-tolerant organisms (Bickler and Buck, 2007; Chapman et al., 2002, 1995; Mandic et al., 2009b; Pollock et al., 2007; Richards, 2011; Wu, 2002).

Several metrics are used as indices of hypoxia tolerance in aquatic organisms. Many species maintain resting rates of oxygen consumption (Ṁ_O2) across a range of high O₂ tensions (P_O2) but will exhibit progressive declines in Ṁ_O2 at low P_O2. The transition from oxyregulation to oxyconformation of Ṁ_O2, which occurs at a P_O2 that is termed the critical O₂ tension (P_crit), is one common measure of hypoxia tolerance (Regan et al., 2019; Richards, 2009; Rogers et al., 2016). P_crit is often calculated as the breakpoint in a two-segmented linear regression (Yeager and Ultsch, 1989). However, studies have found that some species do not exhibit true oxyregulation; some species appear to oxyconform, in which Ṁ_O2 declines progressively as P_O2 falls from normoxic levels (Urbina et al., 2012; Wood, 2018), and the patterns of Ṁ_O2 variation cannot be adequately fitted using two-segmented regression. For this and other reasons, the value of P_crit as an index of hypoxia tolerance has been debated recently (Regan et al., 2019; Wood, 2018). Alternative metrics that summarize the P_O2 dependence of Ṁ_O2, such as the regulation index (a measure of the relative degree of oxyregulation), have been proposed to overcome these criticisms (Mueller and Seymour, 2011; Wood, 2018), but these metrics may represent different physiological information (Regan et al., 2019). Hypoxia tolerance is also reflected by the ability to resist the loss of equilibrium (LOE) in severe hypoxia, as either the time to LOE at a constant level of severe hypoxia (t_LOE) or the P_O2 at LOE (P_LOE) during progressive hypoxia (Borowiec et al., 2016; Crans et al., 2015; Dhillon et al., 2013; Mandic et al., 2013; McBryan et al., 2016). These various indices of hypoxia tolerance have sometimes (Dan et al., 2014; Mandic et al., 2013; Yang et al., 2013), but often not (Crans et al., 2015; Dhillon et al., 2013; Fu et al., 2014; Mathers et al., 2014; Speers-Roesch et al., 2013), been found to co-vary across species.

Fish can respond to and cope with hypoxia by enacting cardiorespiratory adjustments that help maintain cellular O₂ supply, and thus maintain aerobic metabolism, or by reducing cellular O₂ demands via anaerobic metabolism or metabolic depression (Farrell and Richards, 2009; Richards, 2009, 2010). The ability to maintain cellular O₂ supply can obviate the need to reduce cellular O₂ demands, and the relative ability of fish to safeguard cellular O₂ levels should affect their ability to survive in hypoxia. For this reason, indices of hypoxia tolerance that reflect the ability to maintain cellular O₂ supply, such as P_crit, are expected to correlate with those that reflect imminent death, such as t_LOE or P_LOE (McBryan et al., 2016; Speers-Roesch et al., 2013). However, survival in hypoxia could also vary between species depending on their relative ability to reduce cellular O₂ demands, even between species encountering the same reduction in cellular O₂ levels, and in such cases there might be a poor correlation between P_crit and t_LOE or P_LOE (Barnes et al., 2011; Crans et al., 2015; Dhillon et al., 2013; Speers-Roesch et al., 2013). By this rationale, a poor correlation between P_crit and t_LOE or P_LOE could suggest that a reduction in cellular O₂ demands, or the ability to cope with the secondary consequences of such a reduction (e.g. metabolic acidosis), plays a greater role in determining variation in the ability to survive hypoxia exposure.

Variation in hypoxia tolerance results from developmental plasticity, adult phenotypic plasticity and/or evolutionary innovation, in association with variation in many underlying physiological traits. Adult fish can show substantial plasticity that improves hypoxia tolerance in response to chronic hypoxia exposure, in association with morphological and physiological changes in a number of underlying traits involved in oxygen uptake, transport and utilization (Borowiec et al., 2015; Burnett et al., 2007; Du et al., 2016; Fu et al., 2011; Greaney et al., 1980; Martinez et al., 2006; Søllid et al., 2003; Sollid and Nilsson, 2006). Exposure to hypoxia during early development can also can elicit changes in hypoxia tolerance and in other phenotypic traits that persist into adulthood (Blank and Burggren, 2014; Heinrich et al., 2011; Robertson et al., 2014). Enhanced hypoxia tolerance can also evolve over generations, such as when species evolve to become more specialized for life in hypoxia-prone environments, and can be associated with evolved changes in the underlying determinants of O₂ transport and utilization (Hopkins and Powell, 2001; Mandic et al., 2009b; Regan et al., 2017b; Richards, 2011). However, we know much less about the interactions between these processes – namely, the extent to which the plastic responses to chronic hypoxia might differ between species. Some evidence from African cichlids suggests that adaptation to hypoxia-prone swamps can attenuate the changes in brain size in response to developmental hypoxia (Chapman et al., 2008; Crispo and Chapman, 2010), but we know relatively little about the magnitude of interspecific variation in the plasticity of hypoxia tolerance across closely related species.

Killifish from the family Fundulidae (Fig. 1) are well suited for evaluating the roles of phenotypic plasticity and evolutionary innovation in the tolerance of challenging environments. Species of this family are widely distributed and occupy habitats spanning a range of dissolved oxygen, salinity, pH, temperature and other environmental factors (Burnett et al., 2007; Whitehead, 2010; Whitehead et al., 2013). Considerable variation in physiological tolerance of environmental stressors occurs across this family, and even closely related sister taxa can have very different tolerances (Griffith, 1974; Nordlie, 2006; Whitehead, 2010; Whitehead et al., 2013), making Fundulidae a particularly useful model for evolutionary physiology (Burnett et al., 2007). Fundulus heteroclitus exhibits significant plasticity to chronic hypoxia (Borowiec et al., 2015, 2018; Du et al., 2016; Martinez et al., 2006) as well as intraspecific variation in hypoxia tolerance (McBryan et al., 2016), but it is unclear how hypoxia tolerance and its plasticity vary across the family. Therefore, our first two objectives were to characterize (i) the interspecific variation in hypoxia tolerance across Fundulidae and (ii) whether there is interspecific variation in the plasticity of hypoxia tolerance in response to chronic hypoxia. Finally, given recent debate surrounding P_crit and the most appropriate indices of hypoxia tolerance in fishes (Regan et al., 2019; Wood, 2018), our third objective was (iii) to examine whether there is a strong correlation or an uncoupling of different indices of hypoxia tolerance across a variety of taxa and patterns of hypoxia exposure.

Wild populations of killifish were collected with minnow traps, dip nets or seine nets. Fundulus confluentus Goode and Bean 1879 and Fundulus heteroclitus (Linnaeus 1766) were collected from Jekyll Island, GA, USA (31.1039°N, 81.4061°W) with minnow traps. Lucania parva (Baird and Girard 1855) and Lucania goodei Jordan 1880 were collected from the Wakulla River, FL, USA (30.1761°N, 84.245°W) using dip nets. Fundulus diaphanus (Lesueur 1817) and Fundulus rathbuni Jordan and Meek 1889 were collected from Lake Opinicon, ON, Canada (44.559°N, −76.328°W) and Chapel Hill, NC, USA (35.9266474N, −79.0318428W), respectively, with a beach seine. Fish were treated with a 24 h cupramine soak (to treat ectoparasites) followed by 2 weeks of praziquantel treatment (to treat internal worm infections) upon arrival at Louisiana State University, during which they were held at a salinity comparable to that at the collection site (i.e. <0.2 ppt or 12 ppt). Fish were then slowly acclimated to common 0.1 ppt conditions over the course of the next 1–2 weeks and were then held at that salinity until any experimentation began (see below). Experiments were carried out at Louisiana State University for all taxa except F. diaphanus, for which experiments were carried out at McMaster University, and care was taken to provide consistent housing conditions (similar densities, aquarium sizes, etc.) and equivalent experimental treatments across all species at both sites.

At Louisiana State University, fish were initially housed in a filtered recirculating rack system in aerated (normoxic) water (∼20 kPa, 8 mg O₂ l⁻¹) with a salinity of 0.1 ppt, at room temperature (22–26°C) on a 12 h light:12 h dark photoperiod, and were then transferred to 35 l aquaria for 28 day experimental acclimation. At McMaster University, fish were initially housed in 35 l glass aquaria under the conditions described above and were kept in these aquaria for experimental acclimation. Fish at both sites were fed daily with commercial fish pellets (Cargill, Minneapolis, MN, USA). Water chemistry was monitored at least once per week (measurements were similar across tanks and ranged as follows: ammonia <1.0 ppm, nitrates <10 ppm, nitrites <1.0 ppm; pH 8.0–8.3), and water changes were performed as necessary to maintain good water quality. Acclimations to normoxia and hypoxia (see below), respirometry and hypoxia tolerance measurements were conducted in aquaria at least 1 month after arrival at McMaster University or Louisiana State University. All fish were held at room temperature during all acclimations and hypoxia tolerance experiments. Daily temperatures ranged from 22.0 to 26.0°C over the course of the study, but average temperatures for acclimations and hypoxia tolerance experiments did not vary between species or acclimation P_O2. For normoxia and hypoxia acclimations, fish were prevented from accessing the water surface using a plastic grid barrier over which we laid bubble wrap to minimize the diffusion of O₂ from the atmosphere into the water in the aquarium. Because of space limitations, we were restricted to a single tank per treatment, and so did not include tank effects in our analysis. Nevertheless, there were no measurable differences between tanks, with the exception of the desired variation in acclimation P_O2 (hysteresis reported below). All procedures for collecting wild fish and for subsequent experimental treatments were approved by the institutional animal ethics boards of each institution.

A subset of species (F. rathbuni, L. parva, L. goodei) were also exposed for 28 days to constant hypoxia or to nocturnal intermittent hypoxia (12 h hypoxia:12 h normoxia, matched to the photoperiod), as previously described (Borowiec et al., 2015, 2018). Constant hypoxia (P_O2 set point of 2 kPa O₂ with a 0.1 kPa hysteresis) was maintained by bubbling the water with nitrogen gas, mediated by a feedback loop using a fibre-optic oxygen sensor (Loligo Systems, Tjele, Denmark) to regulate the action of a solenoid valve controlling the flow of nitrogen. Intermittent hypoxia was maintained using the same feedback loop as for constant hypoxia, except that hypoxia (2 kPa) was only maintained at night (19:00 h to 07:00 h local time), and the water was bubbled with air to reoxygenate and maintain normoxia during the day. Respirometry and hypoxia tolerance measurements (see below) were conducted after the completion of chronic exposure. Average body mass for fish acclimated to intermittent hypoxia was as follows: F. rathbuni, 1.44±0.22 g and 1.08±0.24 g (t_LOE only); L. parva, 0.41±0.04 g and 0.55±0.08 g (t_LOE only); L. goodei, 0.32±0.02 g and 0.41±0.04 g (t_LOE only). Average body mass for fish acclimated to constant hypoxia was as follows: F. rathbuni, 1.00±0.20 g and 1.26±0.24 g (t_LOE only); L. parva, 0.33±0.05 g and 0.45±0.08 g (t_LOE only); L. goodei, 0.31±0.03 g and 0.34±0.04 g (t_LOE only). See Fig. 4 for the body mass of normoxia-acclimated fish.

Stop-flow respirometry was used to measure resting Ṁ_O2, P_crit, regulation index (RI) and P_LOE (Fig. 2A), using methods consistent with those we have described previously for F. heteroclitus (Borowiec et al., 2015). Fish were habituated overnight in normoxia to a 70 ml cylindrical respirometry chamber that was situated in a large darkened buffer tank. The chamber was connected to a flush pump that circulated water from the buffer tank through the chamber in a ‘flushing circuit’. A second pump circulated water from the chamber in a closed loop (a ‘recirculating circuit’) across a flow-through fibre-optic O₂ sensor (PreSens, Regensburg, Germany). Water flow through both circuits was driven by pumps controlled by AutoResp software (Loligo Systems). Pumps in both the flushing and recirculating circuit were active during the overnight habituation to the respirometry chamber.

Ṁ_O2 measurements began the following morning, and sequential activation and deactivation of the pumps allowed measurement of the change in O₂ concentration due to fish respiration. During flush periods, both pumps were active, the chamber received a steady flow of water from the buffer tank, and no measurements of Ṁ_O2 were conducted. During measurement periods, the flushing pump was deactivated, but continued pumping through the recirculating circuit allowed for measurement of the decline in water oxygen content due to fish respiration. First, resting Ṁ_O2 (indicated by ‘1’ in Fig. 2A) was measured in normoxia. We then measured Ṁ_O2 throughout a progressive stepwise hypoxia protocol, in which the P_O2 of the buffer tank was reduced from ∼20 kPa to 2 kPa in ∼2 kPa steps (∼15 min per step) using the O₂ control system described above, allowing calculation of Ṁ_O2 at 20, 18, 16, 14, 12, 10, 8, 6, 4 and 2 kPa. After Ṁ_O2 was measured at 2 kPa, the chamber was isolated from the buffer tank (by deactivating the flushing circuit) and the fish consumed the remaining O₂ until LOE, and the P_LOE was recorded (indicated by ‘4’ in Fig. 2A). During this period, we recorded Ṁ_O2 at roughly 1.5, 1.0 and 0.5 kPa, but measurements were not possible at each of these nominal P_O2 in every individual (e.g. if the fish had earlier lost equilibrium). A modest number of small fish could not consume O₂ below 0.5 kPa, so we had to open the chamber to the buffer tank and bubble the tank with nitrogen until the fish reached LOE.

Ṁ_O2 was calculated from the change in chamber O₂ concentration over time, as previously recommended (Clark et al., 2013), and is expressed relative to body mass. P_crit was calculated from the Ṁ_O2 and P_O2 data using the R package ‘segmented’ (Muggeo, 2008), which allows the identification and calculation of multiple breakpoints within a single Ṁ_O2–P_O2 curve. To calculate the P_crit of an individual fish, we used the average Ṁ_O2 measurement from 2–3 replicates for the set point P_O2 at each step. The segmented regression model was fitted to a general linear model of the data (Muggeo, 2008). Most individuals exhibited the expected two-segment association between Ṁ_O2 and P_O2, such that P_crit could be calculated as the single breakpoint using the Ṁ_O2 data across all P_O2 (i.e. in the manner described by Yeager and Ultsch, 1989). However, for some individuals from a subset of species, there were two P_O2 breakpoints in the Ṁ_O2–P_O2 relationship (indicated by ‘2’ and ‘3’ in Fig. 2A) rather than a single breakpoint. This pattern resulted from a decline in Ṁ_O2 across a narrow range of P_O2 just below normoxia, followed by a stabilization of Ṁ_O2 across a broader range of intermediate P_O2 (i.e. an absence or appreciable reduction in the slope of decline), and then another phase of more steeply declining Ṁ_O2 (Fig. 2A). There has recently been criticism of the lack of standardized approaches and a tendency for data pruning in calculations of P_crit, possibly in an effort to force a two-segment association on Ṁ_O2–P_O2 data that do not exhibit this pattern of variation (Wood, 2018). With this criticism in mind, we calculated the P_O2 at which each of the two apparent breakpoints occurred, and designated the breakpoint that occurred at the lower P_O2 as the P_crit. For only a very small number of individuals across all treatment groups (N=1 each of F. rathbuni and L. parva), the model could not converge upon a P_crit value as either the single breakpoint or the lower of two breakpoints in the Ṁ_O2–P_O2 relationship, so we do not report P_crit values for these individuals.

We also used the Ṁ_O2 and P_O2 data to calculate RI for each individual, which provides a relative measure of the degree of oxyregulation by comparing the Ṁ_O2 measured across P_O2 with the Ṁ_O2 expected from perfect oxyconformation (Mueller and Seymour, 2011). P_crit and RI are both indices that describe how Ṁ_O2 reacts to changes in P_O2, but they are calculated in different ways and they represent different aspects of this relationship. Whereas P_crit determines the major breakpoint P_O2 in the Ṁ_O2–P_O2 relationship, RI is a metric of the relative degree of oxyconformity across the entirety of the Ṁ_O2–P_O2 relationship from perfect oxyconformity (0) to perfect oxyregulation (1). RI was calculated by first determining the lines of perfect oxyconformation and perfect oxyregulation, which were the lines from the Ṁ_O2 recorded at a P_O2 of ∼20 kPa and the origin (0, 0) or a horizontal line at that Ṁ_O2 to a P_O2 of zero, respectively. RI was then calculated as the area bound by the individual's measured Ṁ_O2–P_O2 relationship and the oxyconformity line (indicated by the shaded region in Fig. 2A), divided by the triangular area bound by the oxyregulation and oxyconformity lines. Therefore, a RI of 1 describes perfect maintenance of the Ṁ_O2 recorded at 20 kPa during hypoxia exposure, whereas a RI of 0 describes perfect oxyconformation.

t_LOE at a sustained level of severe hypoxia was measured in all species (Fig. 2B). Fish were held individually in small chambers (which did not allow access to the surface) within an aquarium of aerated (normoxic) water for an overnight period. During this period, water was continuously circulated through the chambers with aquarium pumps. The following morning, buffer tank P_O2 was rapidly decreased to 0.6 kPa (0.3 mg O₂ l⁻¹), which typically took less than 5–15 min, by the rapid bubbling of nitrogen gas. This P_O2 was held steady until the fish lost equilibrium, and t_LOE was calculated beginning from the time when the P_O2 in the aquaria first reached 0.6 kPa (‘5’ in Fig. 2B). We chose 0.6 kPa because it probably represents a considerable acute hypoxia stressor for even hypoxia-acclimated fish, being well below the P_crit of F. heteroclitus acclimated for 7 days to 2 kPa hypoxia (Borowiec et al., 2015). Accordingly, animals that lost equilibrium before the P_O2 reached 0.6 kPa were assigned a negative t_LOE, which was reflective of the difference in time between when LOE occurred and when the aquarium reached 0.6 kPa.

Data were first checked for normality using a Shapiro–Wilk test (data not shown). For data across normoxia-acclimated fish, we used one-way ANOVA (on ranks in cases when normality was not confirmed) followed by Dunn's or Dunnett’s post hoc tests (as appropriate to the type of ANOVA used), to examine the change in Ṁ_O2 as a function of P_O2 within a species or to compare hypoxia tolerance metrics between taxa. We similarly used two-way ANOVA to test for effects of hypoxia acclimation, species and their interaction. We used least squares linear regression to test for relationships between body mass and Ṁ_O2 or indices of hypoxia tolerance, and for relationships between Ṁ_O2 and indices of hypoxia tolerance. Data from the same set of normoxia-acclimated F. rathbuni, L. parva and L. goodei were used for both the interspecific comparisons (Figs 3–7) and the hypoxia acclimation comparisons (Figs 8 and 9). These statistical analyses were performed using R Studio or GraphPad Prism software (La Jolla, CA, USA).

Phylogenetically independent contrasts were calculated using Mesquite (http://www.mesquiteproject.org) using the PDAP module (http://mesquiteproject.org/pdap_mesquite). We used a previously published, robust maximum likelihood phylogeny of the Fundulidae based on the consensus sequence of the cytochrome b gene for each species (Whitehead, 2010) (Fig. 1). We pruned the tree to only include species for which we measured all character data (species names in bold in Fig. 1), and we used this tree for phylogenetically independent contrasts analysis. Positivized unstandardized contrasts were calculated from absolute data and then standardized by dividing them by the standard deviation of that contrast (e.g. square root of the sum of corrected branch lengths), as previously recommended (Garland et al., 1992; http://www.mesquiteproject.org). We used least squares linear regressions on standardized contrasts to test in a phylogenetically independent manner for the same relationships between body mass, Ṁ_O2 and indices of hypoxia tolerance that are described above.

We quantified the changes in resting Ṁ_O2 during progressive hypoxia (Fig. 3) in 6 species from multiple lineages of fundulid killifish (Fig. 1), in fish that were well acclimated to normoxia in freshwater (0.1 ppt). The response of Ṁ_O2 to declining P_O2 varied between species (Fig. 3), but all species had a significant main effect of P_O2 on Ṁ_O2 (Fig. 3A, F. rathbuni, F_12,171=4.30, P<0.0001; Fig. 3B, F. diaphanus, H₁₃=37.29, P=0.0002; Fig. 3C, L. parva, H₁₃=46.26, P<0.0001; Fig. 3D, L. goodei, H₁₃=33.74, P=0.0007; Fig. 3E, F. heteroclitus, H₁₃=21.61, P=0.0421; Fig. 3F, F. confluentus, H₁₃=34.77, P=0.0005). There was variation in the P_O2 at which Ṁ_O2 first exhibited statistically significant declines relative to Ṁ_O2 in normoxia, with the extremes represented by L. parva, the species that reduced Ṁ_O2 at the highest P_O2 (occurring at ∼10 kPa), and F. heteroclitus and F. confluentus, species that maintained Ṁ_O2 statistically similar to resting Ṁ_O2 until very low (∼0.5 kPa) P_O2.

Along with the variation in the P_O2 at which Ṁ_O2 first decreased from normoxic Ṁ_O2, the general pattern of the Ṁ_O2–P_O2 curve also differed between taxa. For example, in F. diaphanus and L. goodei, Ṁ_O2 changed very little until a low P_O2 was reached, at which point it declined (Fig. 3B,D). Other species, such as F. confluentus and F. heteroclitus, seemed to show a weaker pattern of oxyregulation, with slight declines in Ṁ_O2 at high levels of declining P_O2 (Fig. 3E,F), but nevertheless showed the typical two-segment Ṁ_O2–P_O2 relationship with a single breakpoint P_O2 that we considered to be P_crit. However, some species (e.g. F. rathbuni, L. parva) exhibited an early decline in Ṁ_O2 during mild hypoxia (P_O2 just below normoxia) followed by stabilization of Ṁ_O2 in moderate hypoxia (with the delineation between the two occurring at an upper breakpoint P_O2), and then a second sharper decline in Ṁ_O2 at low P_O2 below P_crit (Fig. 3A,C). In such cases, in which there were three clear segments in the Ṁ_O2–P_O2 curve, we considered the lower breakpoint to represent P_crit.

There was a significant main effect of species on resting Ṁ_O2 in normoxia (H₆=49.64, P<0.0001) and this appeared to be partially driven by the relatively high Ṁ_O2 of L. parva and L. goodei, as well as the low Ṁ_O2 of F. diaphanus compared with those of other taxa (Fig. 3). This variation appeared to result from the relationship between Ṁ_O2 and body mass. As expected, resting Ṁ_O2 in normoxia was negatively correlated with body mass, such that larger animals had lower mass-specific O₂ demands (Fig. 4A), and this relationship remained significant after accounting for evolutionary relatedness between taxa using phylogenetically independent contrasts (Fig. 4B). These relationships remained true even after removing F. diaphanus from the regression in Fig. 4A, or after removing the strong F. diaphanous–F. rathbuni contrast from the regression in Fig. 4B (data not shown).

Indices of hypoxia tolerance varied appreciably between species (Fig. 5). There were highly significant main effects of species on t_LOE (H₆=71.23, P<0.0001; Fig. 5A), P_LOE (H₆=31.34, P<0.0001; Fig. 5B) and P_crit (H₆=26.56, P<0.0001; Fig. 5C), all of which suggested that F. confluentus was the most hypoxia-tolerant species (longest t_LOE, lowest P_LOE, and tied for lowest P_crit), F. rathbuni and F. diaphanus were the least hypoxia-tolerant species (shortest t_LOE, highest P_LOE, highest P_crit), and the other species had intermediate tolerance. RI also varied across species (Fig. 5D), as reflected by a significant main effect of species (H₆=16.98, P=0.0045), but did not follow a similar pattern to P_crit or to the other indices of hypoxia tolerance. RI was generally greater than zero, suggesting that all species exhibited some degree of oxyregulation in hypoxia. Although our different hypoxia tolerance metrics generally agreed with each other (except RI) in species representing the extremes of hypoxia tolerance (i.e. F. confluentus and F. rathbuni), the rank order of species with intermediate levels of tolerance was inconsistent.

Differences in body mass did not appear to make a major contribution to the interspecific variation in hypoxia tolerance. The correlations of t_LOE, P_LOE or P_crit with body mass were not significant for the allometric regressions using either absolute values or phylogenetically independent contrasts (Table 1). The allometric regressions of RI to body mass were significant before (Fig. 6A) and after (Fig. 6C) phylogenetic correction, suggesting that larger species have a higher RI (Table 1). Fish with a lower mass-specific Ṁ_O2 also had a higher RI (Fig. 6B), but this result was only significant after phylogenetic correction (Fig. 6D, Table 1). However, these significant correlations were entirely driven by the single contrast between the sister taxa pair of F. rathbuni and F. diaphanus, which differ appreciably in body mass (Fig. 4A), Ṁ_O2 (Fig. 3A,B) and RI (Fig. 5D).

Given the substantial variation in metrics associated with hypoxia tolerance across Fundulidae (Fig. 5), we examined how different indices of hypoxia tolerance correlated with each other (Table 1). P_LOE and t_LOE exhibited very similar patterns of variation across species (Fig. 5A,B), and the negative correlation between the absolute values of these traits neared significance (P=0.076). There was also a significant positive correlation between the absolute values of P_LOE and P_crit (Fig. 7A), but this correlation was not significant after phylogenetic correction (Fig. 7B, Table 1). RI showed a significant negative correlation with P_crit after accounting for evolutionary relatedness between taxa (Fig. 7D) using phylogenetically independent contrasts, but the correlation of absolute data was not significant (Fig. 7C).

We examined the plasticity of hypoxia tolerance in response to chronic hypoxia exposure in a subset of species – F. rathbuni, L. goodei and L. parva – that appeared to differ in hypoxia tolerance in normoxia (Fig. 5). This was done using two distinct patterns of chronic hypoxia: constant sustained hypoxia and diel cycles of nocturnal hypoxia (‘intermittent hypoxia’, 12 h hypoxia:12 h normoxia, matched closely to the photoperiod). Hypoxia acclimation appeared to affect the qualitative pattern of the response of Ṁ_O2 to declining P_O2, particularly in F. rathbuni and L. parva, such that animals generally reduced Ṁ_O2 less at higher P_O2 (i.e. a less pronounced drop in Ṁ_O2 above the upper P_O2 breakpoint) after acclimation to constant and/or intermittent hypoxia (Fig. 8). The change in the Ṁ_O2–P_O2 curve was most dramatic in F. rathbuni, where the P_O2 at which Ṁ_O2 was first statistically different from Ṁ_O2 at 20 kPa decreased from 8 kPa in normoxia-acclimated fish to ∼0.5 kPa in hypoxia-acclimated fish (Fig. 8A,D,G).

Hypoxia acclimation led to significant improvements in some (though not all) indices of hypoxia tolerance, as reflected by significant main effects of hypoxia acclimation on t_LOE (F_2,117=22.32, P<0.0001), P_LOE (F_2,80=17.73, P<0.0001) and RI (F_2,91=10.10, P=0.0001), but not on P_crit (F_2,89=1.21, P=0.30) (Fig. 9). Although there was some qualitative variation between the responses to constant hypoxia and intermittent hypoxia, both patterns of hypoxia exposure generally appeared to improve most metrics of hypoxia tolerance. These changes were not associated with any significant changes in resting Ṁ_O2, which was unaffected by hypoxia acclimation (F_2,91=1.52, P=0.23), and there was no species×environment interaction for this trait (F_4,91=0.47, P=0.76) (Fig. 9E). As a result, the species differences in resting Ṁ_O2 that were observed in normoxia (Fig. 3) were found to persist across acclimation environments (significant main effect of species, F_2,91=27.52, P<0.0001).

There were species differences in the effects of chronic hypoxia on some indices of tolerance (Fig. 9). The strongest example of these differences was the highly significant species×environment interaction for t_LOE (F_4,117=8.57, P<0.0001), for which there was also a significant main effect of species (F_2,117=95.81, P<0.0001). This significant interaction was driven by the large increase in t_LOE in L. goodei in response to both patterns of chronic hypoxia and the increase in L. parva in response to intermittent hypoxia, with no effects of hypoxia acclimation in F. rathbuni (Fig. 9A). P_crit did not show a significant main effect of species (F_2,89=0.48, P=0.62), but there was a significant species×environment interaction (F_4,89=2.90, P=0.026) driven by the strong reduction in P_crit in F. rathbuni after hypoxia acclimation (Fig. 9C). Neither P_LOE (F_4,80=1.90, P=0.12) nor RI (F_4,91=1.84, P=0.13) had a significant species×hypoxia acclimation interaction, but in both cases, the significant main effects of hypoxia acclimation on hypoxia tolerance appeared to be driven by responses in only one or two of the three species. The effects of hypoxia acclimation on P_LOE were driven by the much stronger reductions in F. rathbuni than in either Lucania species (Fig. 9B). The effects of hypoxia acclimation on RI were driven by increases in both F. rathbuni and L. parva after constant hypoxia, and in only F. rathbuni after intermittent hypoxia (Fig. 9D). These patterns of variation suggest that the responses of each index of hypoxia tolerance to hypoxia acclimation are uncoupled, and there is interspecific variation in the magnitude of these responses.

We compared data across closely related taxa and in response to hypoxia acclimation to investigate the interspecific variation and plasticity of hypoxia tolerance in fundulid killifishes. Ṁ_O2 and various metrics of hypoxia tolerance (t_LOE, P_LOE, P_crit and RI) varied substantially across species (Figs 3 and 5). No single metric could fully describe the variation across species (Fig. 5), and the interspecific variation in different indices of hypoxia tolerance often did not correlate (Fig. 7, Table 1). Hypoxia acclimation generally improved hypoxia tolerance based on some indices (t_LOE, P_LOE and RI), but there was interspecific variation in the magnitude of this plasticity and in the indices of tolerance that were altered by chronic hypoxia (Figs 8 and 9). Our results suggest that hypoxia tolerance is a complex trait that is best appreciated by fully characterizing the Ṁ_O2–P_O2 relationship and by considering multiple indices of hypoxia tolerance.

There was appreciable interspecific variation in hypoxia tolerance (Fig. 5). Relative to that of other teleost fish, the P_crit of fundulid species ranged from comparable to (F. heteroclitus, F. rathbuni, F. diaphanus, L. parva) to far below (F. confluentus, L. goodei) the typical P_crit across various species of fish held at a similar temperature, based on a recent meta-analysis (Rogers et al., 2016). The P_LOE and t_LOE reported here in Fundulus and Lucania are lower than previous measurements in centrarchids (Borowiec et al., 2016; Crans et al., 2015), and are comparable to some previous measurements across common carp and other cyprinids (Dhillon et al., 2013; Fu et al., 2014), suggesting that Fundulidae killifish are relatively tolerant of extreme hypoxia. Among the species examined in this study, F. confluentus was the most hypoxia tolerant overall and the closely related pair of F. rathbuni and F. diaphanus had the weakest tolerance (according to the pattern of variation in t_LOE, P_LOE and P_crit). This may reflect differences in the native habitat between these species: F. rathbuni and F. diaphanus are typically constrained to freshwater streams or lakes (which are often well oxygenated but can become hypoxic, such as during seasonal ice-cover for F. diaphanus), whereas F. confluentus are widely distributed across the highly variable, and frequently hypoxic (Diaz, 2001), estuaries of the southern Atlantic and the Gulf of Mexico (Griffith, 1974; Nordlie, 2006; Whitehead, 2010). Some of the interspecific variation in P_crit could also arise from differences in O₂ supply capacity as a result of interspecific variation in active metabolism and aerobic scope (Farrell and Richards, 2009; Zhang et al., 2018). Better characterization of the habitat occupied by different fundulid species is needed to more clearly understand the relationship between species distribution and hypoxia tolerance in this family, but similar associations have been observed in other taxa. For example, variation in traits associated with hypoxia tolerance are related to differences in distribution between North American Lepomis sunfish, based on the exclusion of the less tolerant bluegill but not the more tolerant pumpkinseed from northern lakes that experience winterkill events (Farwell et al., 2007).

In general, the different indices of hypoxia tolerance often did not correlate with each other (Table 1), as previously observed in other taxa (Crans et al., 2015; Dhillon et al., 2013; Speers-Roesch et al., 2013), probably because these different metrics are partially underlain by distinct physiological mechanisms (Borowiec et al., 2016; Nilsson and Renshaw, 2004). Fish that maintain resting Ṁ_O2 into deeper levels of hypoxia (perhaps from having a higher O₂ transport capacity) would be expected to have lower P_crit and higher RI (Perry et al., 2009; Richards, 2009). If the mechanisms supporting this ability are associated with increases in cellular O₂ levels, they may also contribute to helping fish resist losing equilibrium in hypoxia, which could help explain the correlation between P_crit and P_LOE observed here (Fig. 7A). However, the correlation between the phylogenetically independent contrasts of these traits was not significant, nor were many other correlations between P_LOE or t_LOE and traits expected to have a respiratory component (P_crit and RI) (Table 1). Therefore, the interspecific variation in P_LOE and t_LOE reported here (Fig. 5) is probably explained at least partly by physiological factors independent of the ability to maintain aerobic metabolism in hypoxia.

Variation in the ability to maintain equilibrium in hypoxia could have instead resulted from variation in the relative use of metabolic depression or anaerobic metabolism to help match cellular ATP supply and demand, or in the relative development and sensitivity to metabolic acidosis, rather than variation in the ability to maintain cellular O₂ supply. Metabolic depression can reduce O₂/ATP demands and stretch limited fuel stores through a general reduction in energetically expensive processes in the cell (Guppy and Withers, 1999; Hochachka et al., 1996; Nilsson and Renshaw, 2004). Interspecific variation in the use of metabolic depression could have contributed to some of the observed variation in hypoxia tolerance, although resting Ṁ_O2 in normoxia was not generally observed to be associated with variation in hypoxia tolerance (e.g. the hypoxia-intolerant F. diaphanus had the lowest resting Ṁ_O2 of all species, and L. goodei and L. parva, which had by far the highest resting Ṁ_O2, were intermediate in their hypoxia tolerance). Variation in the capacity and availability of fuel reserves to support anaerobic metabolism could also have been important; in sculpins, for example, interspecific variation in t_LOE is associated with variation in glycogen reserves and lactate dehydrogenase activity in the brain (Mandic et al., 2013; Speers-Roesch et al., 2013). The detrimental effects of metabolic acidosis may have differed between species as well; among triplefin fish, for example, hypoxia-tolerant species are less susceptible to acidosis-induced mitochondrial dysfunction (Devaux et al., 2019).

There was not always a consistent association between P_LOE and t_LOE (Table 1), largely because there appeared to be discordance between these traits in some species (e.g. F. diaphanus had high P_LOE but intermediate t_LOE). This may be reflective of differences in the underlying physiological mechanism of what causes LOE in each situation, at least in some species. For example, t_LOE has been associated with resistance of brain ATP depletion and anaerobic capacity, whereas the mechanisms underlying P_LOE are less well understood (Mandic et al., 2013; Speers-Roesch et al., 2013). A potential cause of such differences may lie in disparities in the rate of hypoxia induction between protocols, which was ∼2 h during the stepwise reductions in P_O2 that led to our measurement of P_LOE but was <15 min for measurement of t_LOE. Variation in the rate of hypoxia induction over the same order of magnitude has been shown to affect P_crit in goldfish, such that slower rates of induction (8 h) were beneficial because they provided more time to reduce P_crit, via increases in gill surface area and haemoglobin–O₂ binding affinity, whereas faster rates of induction (∼0.5–1.5 h) did not (Regan and Richards, 2017). Alternatively, some species may have suffered detrimental effects of the prolonged exposure to moderate hypoxia that occurred in the stepwise hypoxia protocol that was used to determine P_LOE, such that the slower rate of hypoxia induction was a disadvantage for some species compared with others. This last possibility may explain one of the findings reported here, that F. diaphanus were able to maintain equilibrium for nearly 20 min on average in the t_LOE experiment, when P_O2 was rapidly reduced to 0.6 kPa (Fig. 5A), but they exhibited a P_LOE of ∼0.8 kPa in the stepwise hypoxia experiment when hypoxia was induced more slowly. This distinction could have arisen if anaerobic by-products had accumulated over the long period of exposure to stepwise hypoxia, rising to levels that elicited LOE before the fish reached 0.6 kPa O₂.

There was a negative correlation between P_crit and RI, but only after phylogenetic correction (Fig. 7D). The RI has been proposed by some to be a preferable alternative to P_crit (Mueller and Seymour, 2011; Wood, 2018), and we had expected a stronger association between P_crit and RI, given that both of these metrics aim to describe aspects of how hypoxia affects Ṁ_O2. However, while P_crit describes the breakpoint between oxyconformation and oxyregulation, it is possible for Ṁ_O2 to vary above P_crit, which could affect RI with little effect on P_crit. RI and P_crit could be well correlated among species that exhibit only a single breakpoint in the Ṁ_O2–P_O2 relationship. However, the effects of modest hypoxia on Ṁ_O2 that led to the two-breakpoint pattern in the Ṁ_O2–P_O2 relationship in some species may have disrupted the correlation between these metrics. These issues emphasize the value in reporting the full Ṁ_O2–P_O2 relationship to consider how its shape may help explain the relationship (or lack thereof) between RI and P_crit.

RI was the only metric that did not distinguish the most tolerant killifish species (F. confluentus) from the least tolerant (F. rathbuni) (Fig. 5). Our data also suggest that the RI does not adequately describe the Ṁ_O2 responses of some species to hypoxia. For example, the very low RI (∼0.25) of L. parva could be mistakenly interpreted as evidence for oxyconformity of this species, but the Ṁ_O2–P_O2 relationship clearly shows that this species maintains a stable Ṁ_O2 across a broad range of P_O2 from 10 kPa down to its P_crit of ∼3 kPa (Fig. 3). Fundulus rathbuni and F. diaphanus, two closely related species that have similar P_crit and similarly poor hypoxia tolerance (as also reflected by P_LOE and t_LOE), have very different RI values as a result of the ‘upper breakpoint’ pattern seen in F. rathbuni but not in F. diaphanus. Our data suggest that neither P_crit nor RI can adequately represent the nuances of the Ṁ_O2–P_O2 relationship that we observed in some species of killifish, and that the pattern of interspecific variation in RI is largely inconsistent with multiple other indices of hypoxia tolerance. As such, RI may have limited value as a metric of hypoxia tolerance for Fundulus species.

Although some of the distinctions between P_crit and RI probably arise from the upper breakpoint pattern seen in the Ṁ_O2–P_O2 relationship for some species, the cause of this pattern is not entirely clear. One possibility is that some species are especially sensitive to the very subtle effects of disturbance at the beginning of an experiment (e.g. inactivation of the flush pump, etc.), such that resting normoxic Ṁ_O2 was elevated and fish only returned to a stable standard metabolic rate after they became accustomed to these minor disturbances after a few flush-measurement cycles. The species that exhibited this upper breakpoint pattern (L. parva, L. goodei and F. rathbuni) also appeared to be relatively skittish in laboratory conditions and may be more sensitive to minor disturbance, providing some anecdotal support for this possibility. Another possibility for the upper breakpoint pattern in some species is that animals were sensing and responding to minor changes in P_O2 with modest facultative reductions in resting Ṁ_O2 above P_crit. Regardless of the underlying cause of the upper breakpoint pattern, it tends to reduce RI and thus exaggerates the apparent level of oxyconformity, probably without affecting P_crit (if calculated using three-segment regression as done here), P_LOE or t_LOE. Regardless, the R-script used in this study was suitable for modelling both patterns of respirometry data even though the cause of each pattern is unclear.

While this investigation focused on physiological indices of hypoxia tolerance, fish also make important behavioural responses as well, such as the threshold P_O2 and/or rate at which aquatic surface respiration (ASR) occurs. ASR increases as P_O2 declines and greatly enhances survival in hypoxic water (Kramer and McClure, 1982). Interestingly, hypoxia-tolerant sculpins perform ASR at a higher P_O2 than do hypoxia-intolerant sculpins, suggesting ASR may be an important behavioural strategy to survive hypoxia in fishes (Mandic et al., 2009a). While we did not investigate the use of ASR in this study (our fish were prevented from accessing the surface throughout hypoxia acclimation and respirometry experiments), the use of ASR could modulate the P_O2 experienced by fish and thus allow for persistence in hypoxic waters, particularly in species like F. rathbuni that have relatively low hypoxia tolerance. As both behavioural and physiological responses are important for coping with environmental stress in wild fishes, correlating behavioural responses to hypoxia with physiological indices of hypoxia tolerance across fundulids would be a useful area of future study.

Overall, our findings suggest that multiple indices of hypoxia tolerance are needed to appreciate interspecific variation in how fish cope with hypoxia, even when comparing very closely related species, as was done in this study. We agree with recent suggestions (Regan et al., 2019; Wood, 2018) that full characterization of the Ṁ_O2–P_O2 relationship can provide valuable insight that is not represented by the metrics that are calculated from this relationship (i.e. P_crit and RI). P_LOE and t_LOE, indices that are measured directly and reflect the ability to survive hypoxia, are also critical indices of hypoxia tolerance that are often not correlated with P_crit or RI (Table 1) (Speers-Roesch et al., 2013). The critical need to consider multiple metrics comes from evidence that evolutionary variation in the ability to live in hypoxic environments can result from variation in some but not all metrics of hypoxia tolerance. For example, in the example of North American Lepomis sunfish that is discussed above, evolutionary variation in hypoxia tolerance that underlies differences in species distribution in the wild can be explained by interspecific differences in t_LOE but not P_crit (Borowiec et al., 2016; Farwell et al., 2007; Mathers et al., 2014). Among some other species, species differences in distribution can be explained by differences in P_crit (Chapman et al., 2002; Richards, 2011). Furthermore, considering multiple metrics of hypoxia tolerance may provide a more nuanced understanding of the ecological implications of hypoxia. For example, heavy use of anaerobic metabolism below P_crit may prolong survival in hypoxia, but deplete energy reserves for the performance of ecologically relevant traits (e.g. foraging, reproduction, etc.). Given that no single metric can fully explain how fish respond to and cope with hypoxia, and that no single metric can reliably predict interspecific variation in hypoxic niche, we believe that metrics for describing both aerobic respiration and survival in hypoxia should be included in future studies aiming to characterize hypoxia tolerance.

Hypoxia acclimation generally improved hypoxia tolerance, as reflected by the various indices of tolerance measured here, including reduced P_LOE, increased t_LOE and increased RI. Broadening of the P_O2 range for sustaining resting metabolism and/or body posture are common responses to hypoxia acclimation in fishes (Borowiec et al., 2015, 2018; Fu et al., 2011; Regan et al., 2017b; Richards, 2009), and could reflect the combined impact of physiological changes that increase branchial O₂ uptake or circulatory O₂ transport (Matey et al., 2008; Perry et al., 2009; Wells, 2009), adjust the use of anaerobic metabolism (Richards, 2009; Richards et al., 2008; Vornanen et al., 2009) and actively reduce cellular ATP demands (Boutilier, 2001; Hochachka et al., 1996; Richards, 2010). In F. heteroclitus, the mechanisms used to improve hypoxia tolerance after acclimation appear to differ between patterns of hypoxia exposure. Acclimation to constant hypoxia induces a pronounced ∼50% reduction in whole-animal Ṁ_O2, well beyond the reduction observed in naive fish exposed to acute hypoxia, in association with reductions in gill-filament length and in the oxidative capacity of the muscle (Borowiec et al., 2015, 2018). Whether similar mechanisms underlie the improvements in hypoxia tolerance observed here is unclear, insofar as resting Ṁ_O2 was unaffected by acclimation across species (Fig. 9). If these species are capable of metabolic depression, its use may be reserved for more severe levels of hypoxia, as observed in goldfish (Regan et al., 2017a). Acclimation of F. heteroclitus to intermittent hypoxia, by contrast, leads to increases in Ṁ_O2 in hypoxia, and to increases in the oxidative and gluconeogenic capacities of the liver that could hasten the speed of recovery from anaerobic metabolism (Borowiec et al., 2015, 2018). These distinct coping mechanisms enacted by different patterns of hypoxia exposure are both effective at maintaining cellular ATP content and avoiding metabolic acidosis during hypoxia (Borowiec et al., 2018), and improving hypoxia tolerance (Borowiec et al., 2015).

There were interspecific differences in the magnitude of the response to chronic hypoxia, suggesting that there is evolutionary variation in the phenotypic plasticity associated with hypoxia acclimation across fundulid killifishes (Figs 8 and 9). Fundulus rathbuni exhibited greater plasticity than L. parva and L. goodei for several traits (P_LOE, P_crit, RI). One explanation for this may be differences in the signal for plasticity across taxa (e.g. F. rathbuni may experience a greater decrease in tissue P_O2 during hypoxia acclimation). Related to this possibility, the severity of hypoxia relative P_crit could affect the plastic response, as fish were chronically exposed to a P_O2 that was below the normoxic P_crit for two species (F. rathbuni and L. parva) but above it for the other species (L. goodei). This may help explain why P_crit was reduced by exposure to constant hypoxia in F. rathbuni and L. parva but not in L. goodei (Fig. 9C). Another possible explanation is that each species differs in its inherent capacity for plasticity, potentially as a result of interspecific variation in the activation and capacity of different coping mechanisms during exposure to chronic hypoxia. This may explain why L. goodei, which had the lowest P_crit of the three species, showed such a strong improvement in t_LOE in response to constant hypoxia (Fig. 9A). Relatively few detailed investigations into interspecific variation in phenotypic plasticity in response to chronic hypoxia have been done in other fishes, but some previous studies have revealed interspecific variation in how metabolic function and gene expression respond to days of hypoxia exposure (Fu et al., 2014; Mandic et al., 2014). Such interspecific variation in plastic responses to hypoxia probably has important implications for the ecology and distribution of different fish species.

Phenotypic plasticity can facilitate (or sometimes impede) colonization of novel environments, and the magnitude of plasticity can evolve (Fordyce, 2006; Ghalambor et al., 2007; Storz et al., 2010). The observed improvements in hypoxia tolerance associated with hypoxia acclimation probably help fundulid killifish colonize hypoxic environments in the wild. Most of the killifish species studied here were relatively tolerant of hypoxia by standard measures, and many inhabit environments that can become quite hypoxic, including estuaries (F. heteroclitus, F. confluentus, L. parva) and freshwater lakes that experience seasonal ice-cover and hypoxia (F. diaphanus) (Hasler et al., 2009; Nordlie, 2006; Whitehead, 2010). Our findings suggest that hypoxia tolerance is plastic across fundulids, but the manifestation of plasticity can differ between species, which could contribute to the broad distribution of this family across North America (Nordlie, 2006; Whitehead, 2010).

We thank Tiffany Simms Lindsay, Kat Munley, Charlie Brown and other members of the Galvez lab for their assistance. We also thank Catie Ivy for her help with troubleshooting our R script. Matt D. Regan and Jeff G. Richards provided excellent comments on a pre-submission version of the manuscript. We also thank the anonymous reviewers for their comments on a previous version of the manuscript.