Emersion behaviour underlies variation in gill morphology and aquatic respiratory function in the amphibious fish Kryptolebias marmoratus

Oxygen availability in aquatic systems is often spatially and temporally variable (Crispo and Chapman, 2008; Diaz and Rosenberg, 2008). To match metabolic O₂ demands to low O₂ supply in hypoxic conditions, fishes use a variety of behavioural, physiological and morphological adjustments to reduce O₂ requirements and increase O₂ uptake from the water (Richards et al., 2009). For example, fishes may swim away from hypoxic areas or use aquatic surface respiration to ventilate the gills with relatively oxygenated water from the air-water interface (Kramer and McClure, 1982; Chapman and McKenzie, 2009). These behaviours are often complemented by physiological and morphological strategies to increase the capacity for O₂ uptake at the gills (Sollid and Nilsson, 2006; Richards, 2009). Fishes can maximize both functional surface area (the gill area that is perfused with blood), by increasing blood flow, and the total surface area (the maximum gill area that can be perfused), via growth of gill filaments or lamellae (Timmerman and Chapman, 2004; Chapman et al., 2008; Blank and Burggren, 2014). In some species, the total gill surface area can also be increased by reversibly decreasing the size of the interlamellar cell mass (ILCM) – epithelial tissue that fills the water channels between gill lamellae (reviewed in Sollid and Nilsson, 2006; Nilsson, 2007; Nilsson et al., 2012).

Fishes have repeatedly evolved amphibious lifestyles in response to hypoxic aquatic environments (Graham, 1997; Graham and Lee, 2004; Wright and Turko, 2016). Amphibious fishes often respond to severe aquatic hypoxia by leaving the water (emersion) and breathing atmospheric air using extra-branchial air-breathing organs (e.g. Mandic et al., 2009a; Urbina et al., 2011; Regan et al., 2011). In these species, large gill surface areas may thus be unnecessary for maintaining aquatic O₂ uptake, and long gill filaments and lamella may be susceptible to damage or cause increased rates of evaporative water loss during emersion (Nilsson et al., 2012; Wright, 2012). To support the gills during emersion, some amphibious fishes have evolved specialized gill structures (Munshi, 1976). For example, several mudskippers possess stout filaments and lamellae (Low et al., 1988), and the lungfish Protopterus annectens has well-separated lamellae that are reversibly covered with mucus during terrestrial aestivation (Sturla et al., 2002). In the amphibious mangrove rivulus Kryptolebias marmoratus (Poey 1880), gill arches are reversibly stiffened after 7 days of terrestrial acclimation (Turko et al., 2017). Gill remodelling via enlargement of an ILCM also occurs in K. marmoratus out of water (Ong et al., 2007), which may minimize evaporative water loss or provide support for the lamellae, but also reduces total surface area and impairs O₂ uptake when fish return to water (Turko et al., 2012). Thus, total gill surface area in amphibious fishes may be subject to a trade-off in which adequate surface area for aquatic respiration is balanced against the benefits of small gills while out of water.

Our goal was to understand how the environment influences gill remodelling and aquatic respiration of amphibious fishes by comparing wild populations of K. marmoratus from habitats with contrasting O₂ availability. Mangrove rivulus are euryhaline, amphibious killifish that inhabit small ponds and crab burrows throughout mangrove swamps of the tropical western Atlantic (Taylor, 2012; Tatarenkov et al., 2017). These habitats are often hypoxic and high in hydrogen sulphide and, like several other amphibious fishes, mangrove rivulus avoid these stressors by leaving water and breathing through their skin (Taylor, 2012; Wright, 2012; Turko and Wright, 2015). In Belize, we found a population of mangrove rivulus living in a hypoxic (P_O2 ∼3.2 kPa) freshwater pool and other fish inhabiting nearly anoxic (∼0.8 kPa) seawater crab burrows. If the nearly anoxic conditions in the crab burrow habitat causes mangrove rivulus to leave water and breathe air, these fish should have higher rates of emersion, reduced gill surface area (ILCM enlargement) and reduced capacity for aquatic gas exchange. To evaluate these predictions, we measured the critical O₂ tension (P_crit; an indicator of respiratory function), gill morphology and emersion rates in fish from both habitats.

In addition to environmental influences on respiratory function via gill remodelling, intraspecific genetic differences may also influence respiratory performance. Specifically, heterozygous individuals have been hypothesized to be able to mount a more robust response to environmental challenges than inbred, homozygous individuals (a phenomenon known as inbreeding depression; Shull, 1948; Samollow and Soulé, 1983; Shikano and Taniguchi, 2002, Reed and Frankham, 2003). Considering that respiratory function is determined by complex cellular and whole-animal traits involving hundreds of genes, the effects of inbreeding depression may be particularly pronounced (Gracey et al., 2001; Ton et al., 2003; van der Meer et al., 2005; Mandic et al., 2014). However, there are few direct tests of the heterozygote advantage hypothesis in vertebrates, and a meta-analysis of the empirical data suggests that inbreeding effects may be minimal (Chapman et al., 2009). As one of only two known self-fertilizing hermaphroditic vertebrates (Harrington, 1961; Tatarenkov et al., 2009; Avise and Tatarenkov, 2015), K. marmoratus are an excellent model for investigating whether inbreeding reduces respiratory function. In the wild, K. marmoratus hermaphrodites coexist with a smaller number of males, resulting in distinct populations including almost completely homozygous and also largely heterozygous individuals (Mackiewicz et al., 2006a; Turner et al., 2006; Ellison et al., 2011; Tatarenkov et al., 2015). Our second objective was, thus, to determine whether natural variation in heterozygosity (estimated using 32 microsatellite loci) influenced respiratory function (approximated as P_crit) in wild fish.

In December 2012, mangrove rivulus were collected on Long Caye, Lighthouse Reef Atoll, Belize, from three sites within 500 m of each other: a freshwater pond, a cluster of crab burrows and a small pool by the former Calypso resort (Fig. S1). Fish from the freshwater pond and crab burrow sites were compared to determine how the environment influences gill morphology, respiratory function and emersion behaviour, as these two sites differed greatly in abiotic conditions (Table 1). The Calypso site is a designated ‘no take’ site for long-term monitoring of population genetics in mangrove rivulus (R.L.E., D.S.T., A.T., in preparation) and, so, fish from this site were only briefly taken for determination of respiratory function and heterozygosity (fin clips) before being released. Gee minnow traps were used to capture fish in open ponds, and Taylor cup traps were used to sample crab burrows (Davis et al., 2003; Taylor, 1990). Fish were kept in plastic Whirl-Pak bags (60 ml) in water from their respective habitats for up to 48 h before experiments. All experiments were approved by the University of Guelph animal care committee.

Water conditions at each site (temperature, salinity, pH, dissolved O₂) were measured using an electronic field probe (HI 9828, Hanna Instruments, Woonsocket, RI, USA). Aquatic conditions within crab burrows were taken from the same 18 burrows four times over 7 days between 10:00 and 15:00 h. Measurements of the freshwater and Calypso pools were taken at the same location five to six times between 09:00 and 21:00 h, over the span of 8 days. The crab burrow sampling site was located in a dense forest of red mangrove trees (Rhizophora mangle) on the west side of Long Caye, where numerous land crab (Cardisoma guanhumi) burrows were found (Fig. S1). Water in the crab burrows ranged in salinity from brackish to full-strength seawater and was severely hypoxic or anoxic throughout the day (Table 1; Fig. S2). The freshwater pond site contained abundant K. marmoratus and was a shallow (<50 cm depth, ∼200 m²) open pool with plentiful green algae on the northeastern part of Long Caye. The pool contained freshwater (0.3‰) and was moderately hypoxic, with partial pressures of O₂ that ranged from 1.2 kPa in the morning to 4.5 kPa by late afternoon (Table 1; Fig. S2). Among the many K. marmoratus field collection sites described from western Atlantic mangroves, this location is the first known completely freshwater habitat reported for this species (Taylor, 2012). Although it is known that mangrove rivulus can survive in freshwater for several weeks (LeBlanc et al., 2010), to our knowledge, no other field site has exhibited such low salinity.

Aquatic O₂ consumption of size-matched fish (Table 2) was measured using closed glass respirometry chambers (∼50 ml, 28°C), as described previously (Rodela and Wright, 2006; Turko et al., 2012). Briefly, dissolved O₂ (DO) was continually monitored using Clark-type electrodes (Vernier DO-BTA and LabPro, Vernier Software and Technology, Beaverton, OR, USA) connected to a computer running Vernier LoggerPro 3.8 software. Fish were acclimated to the respirometry chambers for 30 min, as preliminary experiments indicated that 30 min was sufficient to allow metabolic rates to stabilize. Chambers were then sealed and O₂ consumption was recorded over a 2–4 h period until DO dropped below 2% of saturation. Background O₂ consumption was measured in chambers without fish for approximately 1 h immediately after each experiment and this value was subtracted from the O₂ consumption values of the fish. Oxygen consumption of freshwater fish (n=11) was measured in well-aerated water obtained from the pond where they were captured. Fish collected from crab burrows (n=10) and the Calypso site (n=11) were tested in relatively clean, well-aerated seawater (34‰). High background rates of O₂ consumption in crab burrow and Calypso pond water due to large amounts of organic matter would have prevented accurate measurements of metabolic rate. Immediately after each trial, fish were killed (using 500 mg l⁻¹ tricaine methanesulphonate buffered to neutral with sodium bicarbonate) and a fin clip was taken for determination of heterozygosity (see below).

Critical O₂ tension (P_crit) was calculated from the O₂ consumption curves using two methods: broken-stick regression (BSR; Yeager and Ultsch, 1989) and nonlinear regression (NLR; Marshall et al., 2013). The traditional BSR approach estimates P_crit as the intersection of the two linear regression lines that best fit the P_O2 versus O₂ consumption plot. BSR estimates of P_crit were calculated using REGRESS, Jeffrey Muday's software (www.wfu.edu/~mudayja/software/o2.exe), which uses the algorithm described in Yeager and Ultsch (1989). The newer NLR technique was designed to calculate P_crit in cases when the transition between O₂ regulation and O₂ conformation is not a clear break-point; instead of using linear regression, this method estimates P_crit from the best nonlinear function that describes the P_O2 versus O₂ consumption relationship (Marshall et al., 2013). Briefly, six separate curves (Michaelis–Menten, power, hyperbola, Pareto and Weibull with or without an intercept) were fitted to a normalized O₂ consumption curve for each fish, and the curve with the lowest corrected Akaike's information criterion (AICc) score was chosen. P_crit was then calculated (slope m=0.065) from the derivative of the chosen function using the equations in Marshall et al. (2013) or using Maple software (Maplesoft, Waterloo, ON, Canada) to solve Weibull functions.

The regulation index (RI) is a complementary measure to P_crit that describes the extent to which an organism regulates O₂ uptake at low environmental P_O2 (Mueller and Seymour, 2011). The RI is calculated from the same O₂ consumption versus P_O2 trace used to calculate P_crit, and represents the relative area under the curve but above the line of equality. Thus, a perfect oxyconformer has an RI value of zero (i.e. O₂ consumption matches P_O2), whereas a hypothetically perfect O₂ regulator maintains O₂ consumption even in near-anoxic conditions (RI=1)_.

Tissue samples were taken immediately from fish collected from both saltwater crab burrows and the freshwater pond. To determine whether the gill morphology of wild-caught fish was phenotypically plastic, additional fish from each of these habitats (n=3–11) were acclimated to a terrestrial environment for 7 or 13 days before being killed. Fish were kept out of water on moist filter paper above a cotton ball reservoir soaked with either seawater (crab burrow fish) or freshwater (freshwater pond fish), as described previously (Ong et al., 2007). Whole heads were fixed in 10% buffered formalin for 24 h, decalcified (Surgipath Decalcifier II, Winnipeg, MB, Canada) for 1 h, and then transferred to 70% ethanol for storage and shipment to Guelph, ON, Canada. Tissues were routinely paraffin embedded, sectioned in 5 µm increments, and stained with haematoxylin and eosin. Photomicrographs were taken using a Nikon Eclipse 90i microscope and measurements were made using ImageJ software (National Institutes of Health, Bethesda, MD, USA). Five gill lamellae per arch were randomly selected, and the length of each lamella and size of the adjacent ILCM was measured (Ong et al., 2007; LeBlanc et al., 2010).

Emersion rates were quantified for individual fish by measuring the proportion of time each individual spent stuck to the side of the holding container above the water (Turko et al., 2011). Fish collected from either crab burrows (n=10) or the freshwater pond (n=10) were placed in translucent 100 ml plastic containers (FisherBrand Collection Containers; Fisher Scientific) containing 60 ml well-aerated water obtained from the collection site. After a 1 h acclimation period, fish were video recorded (Logitech Quickcam Pro, Fremont, CA, USA) for an additional 1 h period. Each container was surrounded by white paper on three sides to prevent fish from observing each other, and care was taken not to disturb the fish for the duration of the acclimation and recording period. The video files were subsequently used to quantify the duration of each emersion event to the nearest minute.

Experimental fish were either killed (freshwater pond and crab burrow sites) or lightly anaesthetized (Calypso collection site) with 500 mg l⁻¹ tricaine methanesulphonate buffered to neutral with sodium bicarbonate, and a small (∼1 mm²) piece of caudal fin was collected using an ethanol-cleaned razorblade. Fin clips were immediately stored in DNA preservative (0.25 mol l⁻¹ ethylenediaminetetraacetic acid, 20% dimethyl sulphoxide, NaCl saturated, pH 7.5) for transport to the University of California, Irvine (Irvine, CA, USA). Besides genotyping fish used in the respiration experiments, we sampled and genotyped fin clips from additional fish in the same locations. The combined samples (i.e. including experimental and non-experimental) were as follows: freshwater pond (n=82), crab burrows site (n=46) and Calypso pond (n=50). These large sample sizes allowed us to evaluate how well experimental fish represented the genetic diversity in our field site. Genetic relatedness and heterozygosity were estimated with 32 nuclear microsatellite loci developed previously for K. marmoratus (Mackiewicz et al., 2006b). DNA preparations, the genotyping protocol and the binning of alleles followed Tatarenkov et al. (2010, 2012). Individual heterozygosity was calculated by summing the number of heterozygous loci in an individual and dividing the sum by the total number of loci. The number of heterozygous loci for each individual was counted using Microsatellite Analyser v. 4.05 (Dieringer and Schlötterer, 2003). Note that the average of individual heterozygosities corresponds to the observed heterozygosity (H_O) of a population. Genetic differences between individuals were estimated with the distance metric based on the proportion of shared alleles (D_PS; Bowcock et al., 1994). Values of D_PS can range from zero (genetically identical) to one (no shared alleles). Calculations of D_PS were done using Microsatellite Analyser. The genetic relationships among individuals were summarized with the neighbour-joining (NJ) trees constructed in PHYLIP software ver. 3.695 (Felsenstein, 1993).

We compared P_crit, RI and emersion behaviour between the crab burrow and freshwater populations of K. marmoratus using Student's t-tests to understand the effect of the environment on these traits. Two-way repeated measures ANOVA and post hoc Holm–Šidák tests were used to compare rates of O₂ consumption between habitats and across environmental P_O2, and two-way ANOVA was used to compare coverage of the gill lamellae between crab burrow and freshwater pond fish directly from the field and after 7 or 13 days of air exposure.

To determine whether homozygosity influenced respiratory function, we used fish collected from all three sites. We first performed simple linear regression analyses with either BSR or NLR estimates of P_crit as the dependent variable and individual heterozygosities as the independent variable. To account for any possible influence of phylogenetic autocorrelation on our results, we also analysed the data using the phylogenetic generalized least squares (pgls) function in the R package caper (https://cran.r-project.org/web/packages/caper/caper.pdf). These analyses were performed using both NJ and unweighted pair group method with arithmetic mean (UPGMA) phylogenies that were built using the 32 microsatellite markers with a Brownian model for evolution at λ values (importance of branch lengths) of 0, 1 or the ‘maximum likelihood’ value calculated by caper.

RStudio (https://www.rstudio.com/) was used to perform the pgls analysis, and SigmaPlot 11 (Systat Software, San Jose, CA, USA) was used for all other analyses (critical α=0.05). Throughout the text, values are given as means±s.e.m.

In normoxic water, the rate of O₂ consumption by mangrove rivulus from the freshwater pond was significantly lower than in fish collected in crab burrows, and this difference persisted at all values of P_O2 above 6 kPa (two-way repeated measures ANOVA interaction, F_1,18=7.633, P<0.001; Fig. 1A). Fish from the freshwater pond had significantly lower P_crit values than those of crab burrow fish, although the magnitude of this effect varied, depending on the technique used to estimate the break-point of the O₂ consumption curve. Using the traditional BSR approach, the P_crit of crab burrow fish was ∼60% higher than that of freshwater pond fish (t₁₉=3.220, P=0.005; Fig. 1B), but using NLR calculations the P_crit of crab burrow fish was only ∼30% higher than that of freshwater pond fish (t₁₉=2.338, P=0.030; Fig. 1B). Estimates of P_crit using BSR and NLR were significantly correlated (linear regression, R²=0.216, F_1,31=8.522, P=0.006) but differed by an average of 2.55±0.30 kPa and by as much as 5.98 kPa (Fig. S3). There was no difference in the RI between fish collected at different locations (t₁₉=1.333, P=0.20; Fig. 1C).

We observed a significant interaction between habitat and acclimation (two-way ANOVA interaction, F_2,41=3.957, P=0.027; Fig. 2) on gill coverage by the ILCM. Fish sampled immediately after capture from crab burrows had significantly more of the interlamellar space covered by an ILCM than fish from the freshwater pond (t=2.786, P=0.008). After 7 days of air exposure, the ILCM had significantly enlarged in both the crab burrow (t=2.476, P=0.018) and freshwater pond fish (t=8.153, P<0.001; Fig. 2C). The ILCM did not enlarge further after an additional 6 days in air in either group (t≤1.22, P>0.2; Fig. 2E). There was no difference in ILCM coverage between freshwater and crab burrow fish after 7 days (t=1.223, P=0.23) or 13 days (t=0.225, P=0.82) of terrestrial acclimation.

Fish collected from crab burrows spent almost 90% of the recording period out of water, significantly more than freshwater pond fish, which emersed ∼30% of the time (t₁₈=3.896, P=0.001; Fig. 3).

Phylogenies using both UPGMA and NJ methods revealed a similar population structure among mangrove rivulus on Long Caye. NJ trees summarizing genetic similarity (and, therefore, relatedness) among fish are shown in Fig. 4A (for experimental animals) and Fig. S4 (for all fish). Both trees illustrate high genetic diversity, evidenced by long branches, although some fish are genetically identical. There were two major lineages of fish present in the freshwater pond (Fig. 4A). Fish within each of these lineages were highly genetically similar. In one freshwater lineage, 37 fish were identical at all 32 loci examined, and another seven fish were distinct at one locus only. There was slightly more genetic diversity in the other freshwater lineage, although they were nonetheless highly similar, being distinct from each other at no more than three of the 32 loci examined. Notably, divergence between the two freshwater lineages (determined by D_PS) is of the same order as the divergence between any random pair of fish captured in crab burrows or the Calypso pool. The crab burrow and Calypso fish effectively form a genetic mosaic with almost no indication of geographic pattern. These fish were far less homozygous (at 77.8±0.05% of microsatellite loci) than freshwater fish (homozygous at 100±0% of microsatellite loci; t₃₀=3.232, P=0.003). Fig. S4 further demonstrates that our experimental fish represent an unbiased sample of the genetic diversity in the study sites; indeed, experimental fish are randomly and evenly scattered among branches of the full genetic tree. Furthermore, the lack of geographic clustering suggests that the physiological patterns described above are not due to idiosyncrasies of particular genetic lineages, but are the result of the environmental conditions in which fish were collected. However, the high degree of relatedness between some individuals within the freshwater pond site suggests that, when working with self-fertilizing hermaphroditic species, it is possible to take a series of replicate measurements on individuals of the same genotype inadvertently. Thus, caution must be exercised when interpreting results.

Heterozygosity did not significantly explain the variation in P_crit calculated using either BSR or NLR, when tested with simple linear regressions (BSR, R²=0.015, F_1,30=0.458, P=0.50; NLR, R²=0.007, F_1,30=0.235, P=0.63; Fig. 4B). Phylogenetically corrected statistical models similarly found no relationship between heterozygosity and either measure of P_crit (all R²≤0.10, all P>0.05; Table 3).

Water-breathing fishes that live in hypoxic habitats are typically characterized by low rates of O₂ consumption and large gill surface areas, which, in turn, reduce P_crit and allow the maintenance of routine aerobic metabolism at low environmental P_O2 (Mandic et al., 2009b). In contrast, we found that wild K. marmoratus captured from nearly anoxic crab burrows had higher rates of O₂ consumption, increased P_crit and reduced total gill surface area relative to fish in a less hypoxic freshwater pool. The probable explanation for these results is that amphibious K. marmoratus emerse far more frequently in poor water conditions, and air exposure is the dominant determinant of gill remodelling and metabolic status in this species. Finally, we found no evidence that P_crit was affected by high or complete homozygosity, indicating the absence of inbreeding depression in these environments.

Gill remodelling is one mechanism used by some fully aquatic fishes to increase total gill surface area and improve O₂ uptake in response to hypoxic environments (Sollid et al., 2003; Dhillon et al., 2013; Tzaneva et al., 2011, 2014). Similarly, mangrove rivulus acclimated to aquatic hypoxia (4.2 kPa) for 7 days in the laboratory had significantly increased gill surface area (smaller ILCMs) compared with normoxic controls, but these fish were forcibly submerged for the duration of the acclimation period (Turko et al., 2012). In wild mangrove rivulus, however, we observed the opposite pattern – fish collected from nearly anoxic crab burrows had significantly enlarged ILCMs (reduced total gill surface area) compared with fish from the freshwater pool. Thus, in wild mangrove rivulus, it appears that the relationship between aquatic P_O2 and gill morphology is mediated by emersion behaviour, such that aquatic hypoxia ultimately results in mangrove rivulus with gills largely covered by an ILCM.

We found that mangrove rivulus captured from crab burrows spent almost 90% of their time out of water, triple the amount of time spent by fish from the freshwater pond. In the laboratory, gill remodelling in mangrove rivulus is caused both by prolonged acclimation (1 week) out of water (Ong et al., 2007; LeBlanc et al., 2010) and by frequent voluntary emersions interspersed with periods submerged in water (Turko et al., 2011). In two isogenic populations of laboratory fish (from Belize and the Bahamas), we previously quantified individual variation in emersion tendencies and found that some mangrove rivulus never left water, whereas others spent 78% of the week-long recording period out of water. The size of the ILCM in these fish was positively correlated with the amount of time each fish spent out of water, but this relationship disappeared after we prevented fish from emersing (with a mesh screen at the air–water interface) for a further week, strongly indicating that emersion behaviour causes enlargement of the ILCM (Turko et al., 2011). Therefore, we think that relatively high rates of emersion in mangrove rivulus from crab burrows versus those from the freshwater pond caused the differences in gill surface area we measured.

High rates of emersion in the crab burrow fish were probably stimulated by the nearly anoxic conditions in this habitat. In the laboratory, acute exposure to aquatic hypoxia induces emersion (Regan et al., 2011; Blewett et al., 2017) and, once the fish are out of water, respiration occurs across the skin (Grizzle and Thiyagarajah, 1987; Cooper et al., 2012) and bucco-opercular cavity (Turko et al., 2014). High concentrations of hydrogen sulphide (Abel et al., 1987) and CO₂ (Robertson et al., 2015) have also been shown to induce emersion in mangrove rivulus. Given the stagnant conditions within the crab burrows, it is plausible that these gasses may have further contributed to the high rates of emersion we observed. In addition to our behavioural assay, we remotely video recorded a crab burrow and observed 11 instances of emersion by mangrove rivulus within 30 min (supplementary movie 1 in Turko and Wright, 2015). Taylor (2012) also reported frequent observations of emersed fish near crab burrows, suggesting that this behaviour is a regular occurrence in the field and not an artefact of the artificial plastic containers used for our behavioural experiments. Anecdotally, fish were never observed emersed at the freshwater pond site when water was present. However, we did not measure emersion behaviour or dissolved O₂ in the freshwater pond at night, when algal respiration may have caused extreme hypoxia and possibly triggered emersion behaviour. If this occurred, presumably the fish did not spend enough time out of water to experience ILCM growth (Turko et al., 2011).

The functional benefits of ILCM enlargement during air exposure are unclear. One hypothesis is that the ILCM supports the lamellae, and prevents collapse and possible fusion of the epithelial tissue (Ong et al., 2007; Nilsson et al., 2012; Wright, 2012). Support of gill arches and filaments in terrestrially acclimated mangrove rivulus is also provided by collagen deposition (Turko et al., 2017), but bony or cartilaginous tissues are not present in lamellae. Proliferation of the ILCM may thus provide an alternative structural mechanism in the absence of buoyant support from water. Alternatively, large ILCMs may reduce evaporative water loss across the gills (Wright, 2012). Considering that rates of bucco-opercular ventilation are low during air exposure (∼5 ventilations per hour; Turko et al., 2014) and mangrove rivulus only survive in highly humid terrestrial habitats, it seems improbable that the primary role of the ILCM is water conservation.

Mangrove rivulus collected from nearly anoxic crab burrows had significantly higher P_crit and rates of O₂ consumption (in normoxia) compared with fish collected from the freshwater pond. These differences are opposite to the typical pattern observed in fully aquatic fishes (e.g. Mandic et al., 2009b), and indicate that amphibious behaviour may shape these aspects of aquatic respiratory function. One possibility is that access to O₂-rich air during frequent emersions enables a larger ‘metabolic engine’ in crab burrow fish (e.g. faster growth, increased reproduction; Biro and Stamps, 2010), which causes the relatively high aquatic metabolic rates we observed. Increased O₂ demand in these fish may also result from physiological costs associated with terrestrial acclimation, such as maintaining enlarged cutaneous ionocytes (LeBlanc et al., 2010). Crab burrow fish were also in significantly worse body condition (Fulton's K) than freshwater pond fish, suggesting that differences in body composition or relative organ masses (e.g. gonads, liver) between the populations may also cause the observed difference in whole-animal metabolic rate (Boldsen et al., 2013). Alternatively, the worse body condition of crab burrow fish may reflect smaller energy reserves caused by the higher metabolic rate. Finally, the relatively high metabolic rate of crab burrow fish may have been caused by larger ionoregulatory demands in seawater than in freshwater. However, acclimation to water of different salinities did not change rates of O₂ consumption in laboratory-reared mangrove rivulus (Turko et al., 2012).

Relatively high rates of normoxic O₂ consumption in crab burrow fish occurred despite reduced gill surface area (enlarged ILCMs) in this population. This is not surprising, as respiratory gas transfer in fish at rest is thought to be perfusion limited and thus does not depend on gill surface area (Perry and Gilmour, 2002, 2010). At rest, only about 60% of the gill lamellae are perfused with blood in rainbow trout, Oncorhynchus mykiss (Booth, 1978). Furthermore, mangrove rivulus can maintain routine rates of O₂ consumption with gill ventilatory frequency as low as five opercular movements per minute (Turko et al., 2012). Finally, O₂ transfer across the skin may also contribute to total O₂ uptake in mangrove rivulus under aquatic conditions (Grizzle and Thiyagarajah, 1987).

In moderately hypoxic water, fishes can maintain O₂ uptake by increasing gill ventilation and perfusing more gill lamellae with blood but, under severe hypoxia, respiratory gas transfer becomes diffusion limited, and large total gill surface area becomes beneficial (Perry and Gilmour, 2002). The enlarged ILCMs of mangrove rivulus from crab burrows may thus impair O₂ uptake under hypoxic aquatic conditions and increase P_crit, as has been demonstrated in Carassius carassius (Sollid et al., 2003) and Fundulus heteroclitus (McBryan et al., 2016). Similarly, enlarged ILCMs in laboratory-reared mangrove rivulus, induced with acclimation to either soft water or terrestrial conditions, were linked to higher P_crit than in control fish with small ILCMs (Turko et al., 2012). In these earlier experiments, there were no differences in normoxic rates of O₂ consumption between treatments. However, in the current study, crab burrow mangrove rivulus had higher normoxic O₂ consumption rates than freshwater pond fish. It is, therefore, possible that the higher values of P_crit we measured in crab burrow fish simply resulted from higher overall O₂ demand, and the enlarged ILCMs of these fish did not meaningfully impair branchial gas exchange. Consistent with this view, the rate of O₂ consumption in crab burrow fish at any environmental P_O2 was never lower than that in freshwater pond fish, despite the difference in ILCM coverage. If enlarged ILCMs impaired branchial gas exchange, absolute rates of O₂ uptake should have been lower in crab burrow fish under hypoxic conditions, all else being equal. One possibility is that compensatory responses help to maintain O₂ uptake, such as increased haemoglobin concentrations and higher O₂ binding affinity of haemoglobin (Turko et al., 2014). Ultimately, both increased O₂ demand and enlarged ILCMs probably contributed to the increased P_crit of crab burrow mangrove rivulus. This was the case in a phylogenetically controlled study of sculpins, where large gill surface area, low metabolic rate and high haemoglobin affinity for O₂ influenced P_crit (Mandic et al., 2009b).

Phenotypically plastic responses to emersion are one probable cause of the differences we observed between fish from crab burrows and those from the freshwater pond, but genetic differences between these sites could also be involved. In support of the plasticity hypothesis, we found that emersion behaviour in wild fish was correlated with ILCM size, and emersion behaviour drove plastic changes in gill morphology as described above (Turko et al., 2011). Furthermore, we found that both freshwater pond and crab burrow fish responded to 1 and 2 weeks of forced air exposure by increasing ILCM coverage to the same final magnitude. These data suggest that the two wild populations had similar scopes for plasticity in ILCM size, consistent with the hypothesis that the variation we observed in wild fish is due to plasticity. Finally, the freshwater pond fish comprised two distinct lineages, and emersion behaviour, metabolic rate, P_crit and gill morphology were similar in the two, as would be expected if phenotypic differences between them were the result of plasticity rather than of genetic divergence. However, it is premature to rule out the hypothesis that genetic differentiation between populations also contributed, as selection could have simply favoured similar phenotypes in each habitat (Stern and Orgogozo, 2009; Losos, 2011). Genetically based phenotypic differences are known to exist between other populations of K. marmoratus (Lin and Dunson, 1995), and cause differences in gill morphology in the closely related F. heteroclitus (McBryan et al., 2016). Common garden and/or reciprocal transplant experiments would be valuable to disentangle environmental and genetic effects.

We were surprised to find that there was no relationship between heterozygosity and P_crit in wild K. marmoratus. As one of only two vertebrate species with a mixed-mating system that results in nearly homozygous (via self-fertilization) or heterozygous (mating with another individual) offspring, mangrove rivulus are an excellent model system for understanding the evolution of sexual reproduction (Harrington, 1961; Tatarenkov et al., 2009; Avise and Tatarenkov, 2015). Typically, heterozygous offspring are assumed to have higher fitness than that of inbred homozygotes (Stearns, 1987). However, guaranteed reproductive success from self-fertilization may outweigh the costs of inbreeding depression, especially when individual animals regularly colonize new habitats (Baker, 1955; Pannell and Barrett, 1998; Avise and Tatarenkov, 2012). We found no evidence of inbreeding depression with respect to P_crit, suggesting that multiple generations of selfing and subsequent selection may have purged deleterious alleles from mangrove rivulus populations (Dolgin et al., 2007; Gimond et al., 2013). Our ability to detect an inbreeding effect may have been restricted by limited statistical power, considering the relatively small sample size (n=32) but, because there was no hint of a relationship between heterozygosity and either measure of P_crit in our data (R²<0.02; P>0.5), we do not think that this is the case. Thus, self-fertilization may allow single mangrove rivulus to colonize unoccupied and diverse habitats within the mangrove forest without paying a respiratory penalty.

Respiratory function in hypoxia is determined by complex phenotypic traits spanning the respiratory cascade (Gracey et al., 2001; Mandic et al., 2014). One might, therefore, expect that heterozygosity provides a greater diversity of alleles that could, in turn, enhance respiratory function, but this was not the case in our study. Consistent with our results, homozygosity did not impair developmental stability in nine populations of mangrove rivulus (Taylor, 2001). However, more heterozygous K. marmoratus have lower parasite loads, suggesting a cost to inbreeding in some circumstances (Ellison et al., 2011). Resistance to parasites is conferred by only a small set of major histocompatibility complex genes (Ellison et al., 2012). In contrast, the complexity of the respiratory system may protect against inbreeding (Wagner, 2005; Joyner, 2013). For example, mangrove rivulus can modify gill surface area (Ong et al., 2007), gill ventilation (Turko et al., 2012), blood flow (Cooper et al., 2012), haemoglobin concentration and binding affinity (Turko et al., 2014), and/or O₂ delivery (Brunt et al., 2016). Modifications at one or several of these steps of the O₂ transport cascade may be able to compensate for deleterious effects of inbreeding (if they are present) at any of the other levels (Richards et al., 2009). Furthermore, even highly inbred fish (as estimated using microsatellites or known from the pedigree of the laboratory lineages) are not completely homozygous when the whole genome is considered. For example, Lins et al. (2018) showed that genomes of inbred K. marmoratus lineages harbour 0.031–0.055% heterozygous sites, which amounts to tens of thousands of single nucleotide polymorphisms. At least some of this heterozygosity may be physiologically important and preserved by natural selection.

Understanding how animals interact with their environment is crucial for explaining patterns of physiological acclimatization and adaptation. In wild mangrove rivulus, we found that animals from nearly anoxic, saline crab burrows had higher rates of O₂ consumption, increased P_crit and less gill surface area than those from a moderately hypoxic freshwater pond, the opposite pattern to that generally observed in fully aquatic fishes. Instead, mangrove rivulus that inhabit crab burrows are probably highly amphibious and use the hypoxic burrow water only as a short-term refuge in response to disturbance. This largely amphibious habit may, in turn, cause gill remodelling and high metabolic rates, which then prevent these fish from maintaining routine rates of O₂ consumption in hypoxic conditions (increased P_crit). Aquatic hypoxia has been hypothesized to be an important factor promoting the evolution of terrestriality both in early tetrapods and in extant amphibious fishes (Graham, 1997). Our results provide further evidence for this hypothesis by suggesting that hypoxic aquatic conditions and emersion behaviour may interact and cause positive feedback that encourages increasingly terrestrial lifestyles.

We thank the staff of Itza Lodge for logistical support, Michael Wells for analysing emersion videos, Joseph Turko for mathematical assistance and two anonymous reviewers for helpful comments.